SLiMv2.1: SLiMv2.1/derivs.f90 Source File

00001 !================================================================
00002 subroutine derivs(PE,KE)
00003   ! Calculates RHS
00004   implicit none     
00005   include 'param.f90'
00006   include 'terms3d.f90'           
00007   include 'functionlist.f90'
00008   complex fdtemp(nx/2+1,ny,nz_max,nvar)
00009   real temprho(nx,ny),tempbx(nx,ny),tempby(nx,ny),dt2(nx,ny,3)
00010   real  norm,PE,KE,PE1,KE1,fac(nx,ny,nz),fac1
00011   integer i,j,k,n,nonstep,flow_flag
00012   parameter(norm=1./float(nx*ny)) 
00013   !FFT IS NOT NORMALIZED  
00014   nonstep=3-step !1,2 -> 2,1
00015   fdtemp=fd
00016 
00017   ! dr -> RS
00018   call sfftw_execute(plan_f2rd_1)
00019   call sfftw_execute(plan_f2rd_2)
00020   call sfftw_execute(plan_f2rd_3)
00021   call sfftw_execute(plan_f2rd_4)
00022   call sfftw_execute(plan_f2rd_5)
00023   call sfftw_execute(plan_f2rd_6)
00024   fd=fdtemp
00025   rd=rd*norm
00026 
00027 
00028   !Diagnostics not real energies
00029   PE1=sum(rd(:,:,3:nz-2,3)*coeff(:,:,3:nz-2,1,2))
00030   KE1=sum((rd(:,:,3:nz-2,4)**2+rd(:,:,3:nz-2,5)**2+rd(:,:,3:nz-2,6)**2)&
00031        &                                          *coeff(:,:,3:nz-2,1,2))
00032   call mpi_reduce (PE1, PE, 1, MPI_REAL, MPI_SUM, 0, MPI_COMM_WORLD, ierr )
00033   call mpi_bcast (PE, 1, MPI_REAL, 0, MPI_COMM_WORLD, ierr)
00034   call mpi_reduce (KE1, KE, 1, MPI_REAL, MPI_SUM, 0, MPI_COMM_WORLD, ierr )
00035   call mpi_bcast (KE, 1, MPI_REAL, 0, MPI_COMM_WORLD, ierr)
00036   !KE=KE1
00037   !PE=PE1
00038 !!!!FIRST FIND THE EULERIAN DENSITY PERTURBATION!!!!!!!!!!!!!!!
00039   !===================FIND RHO' on cell edges========================
00040   ! rho_0 dr :RS
00041   rvtemp(:,:,:,1)=coeff(:,:,:,1,nonstep)*rd(:,:,:,1)
00042   rvtemp(:,:,:,2)=coeff(:,:,:,1,nonstep)*rd(:,:,:,2)
00043   rvtemp(:,:,:,3)=coeff(:,:,:,1,nonstep)*rd(:,:,:,3)
00044   ! d/dz(rho_0 dr_z) :RS 
00045   rvtemp(:,:,:,3)=ddze_div(rvtemp(:,:,:,3))
00046   ! rho_0 dr_x, rho_0 dr_y -> FS
00047   call sfftw_execute(plan_r2fvtemp_1)
00048   call sfftw_execute(plan_r2fvtemp_2)
00049   !  call sfftw_execute(plan_r2fvtemp_3)
00050 !!!!! NOW STORE FD/dens in fdtemp 
00051   !  fdtemp=fvtemp 
00052   ! d/dx(rho_0 dr_x) :FS
00053   fvtemp(:,:,:,1)=ddxe(fvtemp(:,:,:,1))
00054   ! d/dy(rho_0 dr_y) :FS       
00055   fvtemp(:,:,:,2)=ddye(fvtemp(:,:,:,2))
00056   ! d/dx(rho_0 dr_x) ->RS
00057   call sfftw_execute(plan_f2rvtemp_1)
00058   rvtemp(:,:,:,1)=rvtemp(:,:,:,1)*norm
00059   ! d/dy(rho_0 dr_y) ->RS
00060   call sfftw_execute(plan_f2rvtemp_2)
00061   rvtemp(:,:,:,2)=rvtemp(:,:,:,2)*norm
00062   rrhoe=-(rvtemp(:,:,:,1)+rvtemp(:,:,:,2)+rvtemp(:,:,:,3))
00063 
00064   !rrho is evaluated at edge
00065 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00066 !!!!!!NEXT FIND THE EULERIAN PRESSURE PERTURBATION!!!!!!!!!!!!
00067   ! again at edge
00068   !(GAMMA_1 P_0/rho_0) rho_1 + (GAMMA_1 P_0 Grad(rho_0)/rho_0 - Grad(P_0)) dr :RS  
00069   !WORK most of noise is in rvtemp(:,:,:,3)
00070   rPress= coeff(:,:,:,5,nonstep)*rrhoe  &
00071        &   +coeff(:,:,:,6,nonstep)*rd(:,:,:,1)&
00072        &   +coeff(:,:,:,7,nonstep)*rd(:,:,:,2)&
00073        &   +coeff(:,:,:,8,nonstep)*rd(:,:,:,3)
00074   !-------------------------------------------------------------
00075   !-------------------dB-----------------------------------------
00076   !
00077   ! dB=curl{(dx) cross B_0}
00078   ! at centre
00079   !first calculate (dx) cross B_0 in real space
00080   rvtemp(:,:,:,1)= rd(:,:,:,2)*coeff(:,:,:,17,nonstep)&
00081        &          -rd(:,:,:,3)*coeff(:,:,:,16,nonstep)
00082   rvtemp(:,:,:,2)= rd(:,:,:,3)*coeff(:,:,:,15,nonstep)&
00083        &          -rd(:,:,:,1)*coeff(:,:,:,17,nonstep)
00084   rvtemp(:,:,:,3)= rd(:,:,:,1)*coeff(:,:,:,16,nonstep)&
00085        &          -rd(:,:,:,2)*coeff(:,:,:,15,nonstep)
00086 
00087   ! next calculate dB
00088   call sfftw_execute(plan_r2fvtemp_1)
00089   call sfftw_execute(plan_r2fvtemp_2)
00090   call sfftw_execute(plan_r2fvtemp_3)
00091   !------x component-----------------  
00092   fdb(:,:,:,1)=ddyc(fvtemp(:,:,:,3))
00093   call sfftw_execute(plan_f2rdb_1)
00094   rdb(:,:,:,1)= rdb(:,:,:,1)*norm-ddzc_curl(rvtemp(:,:,:,2))
00095   !------y component----------------- 
00096   fdb(:,:,:,2)=-ddxc(fvtemp(:,:,:,3))
00097   call sfftw_execute(plan_f2rdb_2)
00098   rdb(:,:,:,2)= rdb(:,:,:,2)*norm+ddzc_curl(rvtemp(:,:,:,1))
00099   !------z component-----------------  
00100   fdb(:,:,:,3)=ddxc(fvtemp(:,:,:,2))-ddyc(fvtemp(:,:,:,1))
00101   call sfftw_execute(plan_f2rdb_3)
00102   rdb(:,:,:,3)= rdb(:,:,:,3)*norm
00103 
00104 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00105   if (istop) then
00106      !===========================BCs=======================================
00107      !top boundary-Cally-like
00108      !--------Magnetic Field
00109      ! as per Cameron, Gizon and Daiffallah (2007)
00110 !!!!!!!!!!FIRST LEVEL N-1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00111      ! Bz is known, so first find bx, by then P and finally rho
00112      tempbx = rdb(:,:,nz-1,1)
00113      tempby = rdb(:,:,nz-1,2)
00114      temprho= rrhoe(:,:,nz-1)
00115      rdb(:,:,nz-1,1)=-rdb(:,:,nz-1,3)*coeff(:,:,nz-1,37,nonstep) &
00116           & -rd(:,:,nz-1,1)*coeff(:,:,nz-1,42,nonstep)           &
00117           & -rd(:,:,nz-1,2)*coeff(:,:,nz-1,43,nonstep)           &
00118           & -rd(:,:,nz-1,3)*coeff(:,:,nz-1,44,nonstep)
00119      rdb(:,:,nz-1,2)=-rdb(:,:,nz-1,3)*coeff(:,:,nz-1,38,nonstep) &
00120           & -rd(:,:,nz-1,1)*coeff(:,:,nz-1,45,nonstep)           &
00121           & -rd(:,:,nz-1,2)*coeff(:,:,nz-1,46,nonstep)           &
00122           & -rd(:,:,nz-1,3)*coeff(:,:,nz-1,47,nonstep)
00123 
00124      !     rpress(:,:,nz-1)=-( rd(:,:,nz-1,1)*coeff(:,:,nz-1,39,nonstep)           &
00125      !          &             +rd(:,:,nz-1,2)*coeff(:,:,nz-1,40,nonstep)           &
00126      !          &             +rd(:,:,nz-1,3)*coeff(:,:,nz-1,41,nonstep))          &
00127      !          &           -( coeff(:,:,nz-1,15,nonstep)*rdB(:,:,nz-1,1)          &
00128      !          &             +coeff(:,:,nz-1,16,nonstep)*rdB(:,:,nz-1,2)          &
00129      !          &             -coeff(:,:,nz-1,17,nonstep)*rdB(:,:,nz-1,3))/(4*pi) 
00130      !correction for Lorentz reduction
00131      rpress(:,:,nz-1)=-( rd(:,:,nz-1,1)*coeff(:,:,nz-1,39,nonstep)           &
00132           &             +rd(:,:,nz-1,2)*coeff(:,:,nz-1,40,nonstep)           &
00133           &             +rd(:,:,nz-1,3)*coeff(:,:,nz-1,41,nonstep))          &
00134           &           -( coeff(:,:,nz-1,15,nonstep)*rdB(:,:,nz-1,1)          &
00135           &             +coeff(:,:,nz-1,16,nonstep)*rdB(:,:,nz-1,2)          &
00136           &             -coeff(:,:,nz-1,17,nonstep)*rdB(:,:,nz-1,3))/(4*pi)* &
00137           &      coeff(:,:,nz-1,1,nonstep)**2*coeff(:,:,nz-51,12,nonstep)**2  
00138 
00139      ! temprho is new density at (n-1)
00140      rrhoe(:,:,nz-1)=(      rpress(:,:,nz-1)                                 &
00141           &            +rd(:,:,nz-1,1)*coeff(:,:,nz-1,9,nonstep)             & 
00142           &            +rd(:,:,nz-1,2)*coeff(:,:,nz-1,10,nonstep)            &
00143           &            +rd(:,:,nz-1,3)*coeff(:,:,nz-1,11,nonstep))           &
00144           &            /coeff(:,:,nz-1,5,nonstep)                            &
00145           &            -(rd(:,:,nz-1,1)*coeff(:,:,nz-1,2,nonstep)            & 
00146           &             +rd(:,:,nz-1,2)*coeff(:,:,nz-1,3,nonstep)            &
00147           &             +rd(:,:,nz-1,3)*coeff(:,:,nz-1,4,nonstep))  
00148      !can now compute new rd(:,:,nz,1:2)  
00149      rd(:,:,nz,3)=rd(:,:,nz,3)-2.*dz*coeff(:,:,nz-1,12,1)*(rrhoe(:,:,nz-1)-temprho)  
00150      !     rd(:,:,nz,1)=(-2*dz*(tempbx-rdb(:,:,nz-2,1))+coeff(:,:,nz,15,nonstep)*rd(:,:,nz,3)) &
00151      !          &              *coeff(:,:,nz,47,nonstep)
00152      !     rd(:,:,nz,2)=(-2*dz*(tempbx-rdb(:,:,nz-2,2))+coeff(:,:,nz,16,nonstep)*rd(:,:,nz,3)) &
00153      !          &              *coeff(:,:,nz,47,nonstep)
00154 
00155      !     dt2(:,:,1:3)=dt2+rd(:,:,nz,1:3)
00156      !     if(step .eq. 1) dt2=dt2/(dt)
00157      !     if(step .eq. 2) dt2=dt2/(2.*dt)
00158      !     rd(:,:,nz,1:3)=dt2
00159 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00160      !NOTE HAVE NOW DISPLACMENT AT z=nz
00161      ! would like to have the correct P,B,rho there as well.
00162      !for B only need to work out rdb(:,:,nz,3) : then we also have P and hence rho.
00163 
00164      rdb(:,:,nz,1)=-rdb(:,:,nz,3)*coeff(:,:,nz,37,nonstep) &
00165           & -rd(:,:,nz,1)*coeff(:,:,nz,42,nonstep) &
00166           & -rd(:,:,nz,2)*coeff(:,:,nz,43,nonstep) &
00167           & -rd(:,:,nz,3)*coeff(:,:,nz,44,nonstep)
00168      rdb(:,:,nz,2)=-rdb(:,:,nz,3)*coeff(:,:,nz,38,nonstep) &
00169           & -rd(:,:,nz,1)*coeff(:,:,nz,45,nonstep) &
00170           & -rd(:,:,nz,2)*coeff(:,:,nz,46,nonstep) &
00171           & -rd(:,:,nz,3)*coeff(:,:,nz,47,nonstep)
00172      rdb(:,:,nz,:)=0
00173      rpress(:,:,nz)=-( rd(:,:,nz,1)*coeff(:,:,nz,39,nonstep)            & 
00174           &           +rd(:,:,nz,2)*coeff(:,:,nz,40,nonstep)            &
00175           &           +rd(:,:,nz,3)*coeff(:,:,nz,41,nonstep))           & 
00176           &         -( coeff(:,:,nz,15,nonstep)*rdB(:,:,nz,1)           &
00177           &           +coeff(:,:,nz,16,nonstep)*rdB(:,:,nz,2)           &
00178           &           -coeff(:,:,nz,17,nonstep)*rdB(:,:,nz,3))/(4*pi)    
00179      rrhoe(:,:,nz)=(   rpress(:,:,nz)                                    &
00180           &            +rd(:,:,nz,1)*coeff(:,:,nz,9,nonstep)          & 
00181           &            +rd(:,:,nz,2)*coeff(:,:,nz,10,nonstep)         &
00182           &            +rd(:,:,nz,3)*coeff(:,:,nz,11,nonstep))        &
00183           &            /coeff(:,:,nz,36,step)*coeff(:,:,nz,1,nonstep) &
00184           &            -(rd(:,:,nz,1)*coeff(:,:,nz,2,nonstep)         & 
00185           &             +rd(:,:,nz,2)*coeff(:,:,nz,3,nonstep)         &
00186           &             +rd(:,:,nz,3)*coeff(:,:,nz,4,nonstep)) 
00187 
00188   endif
00189 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00190   !bottom boundary
00191   if(isbottom) then
00192      !  rrhoe(:,:,1)=2*rrhoe(:,:,2)-rrhoe(:,:,3)
00193      rdB(:,:,1,:)=rdB(:,:,2,:)
00194      !     !  rdB(:,:,1,:)=0 
00195   endif
00196   !===============================================
00197 
00198 
00199 
00200   !===================FIND RHO' on cell centres=========
00201   ! rho_0 dr :RS
00202   rvtemp(:,:,:,1)=coeff(:,:,:,1,nonstep)*rd(:,:,:,1)
00203   rvtemp(:,:,:,2)=coeff(:,:,:,1,nonstep)*rd(:,:,:,2)
00204   rvtemp(:,:,:,3)=coeff(:,:,:,1,nonstep)*rd(:,:,:,3)
00205   ! d/dz(rho_0 dr_z) :RS 
00206   rvtemp(:,:,:,3)=ddzc_div(rvtemp(:,:,:,3))
00207   ! rho_0 dr_x, rho_0 dr_y -> FS
00208   call sfftw_execute(plan_r2fvtemp_1)
00209   call sfftw_execute(plan_r2fvtemp_2)
00210   !  call sfftw_execute(plan_r2fvtemp_3)
00211 !!!!! NOW STORE FD/dens in fdtemp 
00212   !  fdtemp=fvtemp 
00213   ! d/dx(rho_0 dr_x) :FS
00214   fvtemp(:,:,:,1)=ddxc(fvtemp(:,:,:,1))
00215   ! d/dy(rho_0 dr_y) :FS       
00216   fvtemp(:,:,:,2)=ddyc(fvtemp(:,:,:,2))
00217   ! d/dx(rho_0 dr_x) ->RS
00218   call sfftw_execute(plan_f2rvtemp_1)
00219   rvtemp(:,:,:,1)=rvtemp(:,:,:,1)*norm
00220   ! d/dy(rho_0 dr_y) ->RS
00221   call sfftw_execute(plan_f2rvtemp_2)
00222   rvtemp(:,:,:,2)=rvtemp(:,:,:,2)*norm
00223   rrhoc=-(rvtemp(:,:,:,1)+rvtemp(:,:,:,2)+rvtemp(:,:,:,3))
00224 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00225   !BC
00226   if(istop) then
00227      !  rrhoc(:,:,nz-1)=(rrhoe(:,:,nz-1)+rrhoe(:,:,nz))/2.
00228      rrhoc(:,:,nz)=rrhoc(:,:,nz-1)
00229   endif
00230   !
00231   rpress(:,:,1)=2.*rpress(:,:,3)-rpress(:,:,4)
00232   rpress(:,:,2)=2.*rpress(:,:,3)-rpress(:,:,4)
00233   rrhoe(:,:,1)=2.*rrhoe(:,:,2)-rrhoe(:,:,3)
00234   rrhoc(:,:,1)=2.*rrhoc(:,:,2)-rrhoc(:,:,3)
00235 
00236   !!
00237 
00238   !Copy to Fourier space (only)
00239   call sfftw_execute(plan_r2frhoe)
00240   call sfftw_execute(plan_r2frhoe)
00241   call sfftw_execute(plan_r2fpress)
00242   call sfftw_execute(plan_r2fdb_1)
00243   call sfftw_execute(plan_r2fdb_2)
00244   call sfftw_execute(plan_r2fdb_3)
00245   !==========================================================
00246   !-------------------------------------------
00247   !----------------dJ-------------------------
00248   fdJ(:,:,:,1)=ddye(fdB(:,:,:,3))
00249   call sfftw_execute(plan_f2rdJ_1)
00250   rdJ(:,:,:,1)=rdJ(:,:,:,1)*norm-ddze_curl(rdB(:,:,:,2))
00251   !-----y component---------------
00252   fdJ(:,:,:,2)=-ddxe(fdB(:,:,:,3))
00253   call sfftw_execute(plan_f2rdJ_2)
00254   rdJ(:,:,:,2)=rdJ(:,:,:,2)*norm+ddze_curl(rdB(:,:,:,1))
00255   !-----z component------------- 
00256   fdJ(:,:,:,3)=ddxe(fdB(:,:,:,2))-ddye(fdB(:,:,:,1))
00257   call sfftw_execute(plan_f2rdJ_3)
00258   rdJ(:,:,:,3)=rdJ(:,:,:,3)*norm
00259 
00260   !-----------------------------------------------
00261   !-----------Next the Lorentz force--------------
00262   !first calculate (dJ) cross B_0 in real space
00263   rLor(:,:,:,1)= rdJ(:,:,:,2)*coeff(:,:,:,17,step)&
00264        &        -rdJ(:,:,:,3)*coeff(:,:,:,16,step)
00265   rLor(:,:,:,2)= rdJ(:,:,:,3)*coeff(:,:,:,15,step)&
00266        &        -rdJ(:,:,:,1)*coeff(:,:,:,17,step)
00267   rLor(:,:,:,3)= rdJ(:,:,:,1)*coeff(:,:,:,16,step)&
00268        &        -rdJ(:,:,:,2)*coeff(:,:,:,15,step)
00269   !next add calculate (J_0) cross dB in real space
00270   rLor(:,:,:,1)=rLor(:,:,:,1)+&
00271        & coeff(:,:,:,19,step)*rdb(:,:,:,3)-coeff(:,:,:,20,step)*rdb(:,:,:,2)
00272   rLor(:,:,:,2)=rLor(:,:,:,2)+&
00273        &coeff(:,:,:,20,step)*rdb(:,:,:,1)-coeff(:,:,:,18,step)*rdb(:,:,:,3)
00274   rLor(:,:,:,3)=rLor(:,:,:,3)+&
00275        &coeff(:,:,:,18,step)*rdb(:,:,:,2)-coeff(:,:,:,19,step)*rdb(:,:,:,1)
00276 
00277 
00278   !CGS normalization
00279   rlor=rlor/(4.*3.1415926535)
00280   !Scale Lorentz force above photosphere
00281   !fac=exp((nz_ghostglobal-80-kk(k))/8.)
00282   !Va/Cs
00283   fac=(coeff(:,:,:,15,step)**2+coeff(:,:,:,16,step)**2+coeff(:,:,:,17,step)**2)/(4.*pi) &
00284        &             *coeff(:,:,:,12,step)/coeff(:,:,:,5,step)
00285   fac=1-(fac/5.)/(1.+(fac/5.))
00286 
00287   !next scales B so that Va will scale up to a maimum of 5*Cs
00288   !With next to no change around Cs=Va
00289   do i=1,3
00290      do k=1,nz
00291         if(kk(k) .ge. nz_ghostglobal-120) rlor(:,:,k,i)=rlor(:,:,k,i)*fac(:,:,k) ! &
00292         !*coeff(1,1,k,1,step)/1.9e-7
00293      enddo
00294   enddo
00295   do i=1,3
00296      do k=1,50
00297         fac1=exp(-(51.-kk(k))/10.)
00298         if(kk(k) .le. 50) rlor(:,:,k,i)=rlor(:,:,k,i)*fac1
00299      enddo
00300   enddo
00301 
00302   ! kinematic
00303   !  rLor=0
00304 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00305 !!!!!!!!!!!!!!!!!!!!!!!!!!!!NO FLOW!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00306 
00307 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00308 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00309   flow_flag=0
00310   if(flow_flag .eq.1) then
00311   write(*,*) "WARNING ASSUMES CARTESIAN"
00312      !---------------------------------------------------------------------------
00313      !------------------------Calculate dU---------------------------------------
00314      ! U'= d/dt dr + (U_0 . nabla ) dr - (dr . nabla) U_0
00315 
00316      ! first centred
00317      rdUc(:,:,:,1)=centre(rd(:,:,:,1))
00318      rdUc(:,:,:,2)=centre(rd(:,:,:,2))
00319      rdUc(:,:,:,3)=centre(rd(:,:,:,3))
00320      !first x_component
00321      ! + (U_0 . nabla ) dr
00322      !calculate d/dx_i rd_j  and move into real space
00323      fvtemp(:,:,:,1)=ddxc(fd(:,:,:,1))
00324      fvtemp(:,:,:,2)=ddyc(fd(:,:,:,1))
00325      call sfftw_execute(plan_f2rvtemp_1)
00326      call sfftw_execute(plan_f2rvtemp_2)
00327      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00328      rvtemp(:,:,:,3)=ddzc(rd(:,:,:,1))
00329      rduc(:,:,:,1)=rduc(:,:,:,1)+coeff(:,:,:,21,step)*rvtemp(:,:,:,1) &
00330           &                     +coeff(:,:,:,22,step)*rvtemp(:,:,:,2) &
00331           &                     +coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00332      ! - (dr . nabla ) U_0
00333      rduc(:,:,:,1)=rduc(:,:,:,1)-centre(rd(:,:,:,1))*coeff(:,:,:,24,nonstep)&
00334           &                     -centre(rd(:,:,:,2))*coeff(:,:,:,27,nonstep)&
00335           &                     -centre(rd(:,:,:,3))*coeff(:,:,:,30,nonstep)
00336      !end x_component
00337      ! second y_component
00338      ! + (U_0 . nabla ) dr
00339      !calculate d/dx_i rd_j  and move into real space
00340      fvtemp(:,:,:,1)=ddxe(fd(:,:,:,2))
00341      fvtemp(:,:,:,2)=ddye(fd(:,:,:,2))
00342      call sfftw_execute(plan_f2rvtemp_1)
00343      call sfftw_execute(plan_f2rvtemp_2)
00344      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00345      rvtemp(:,:,:,3)=ddze(rd(:,:,:,2))
00346      rduc(:,:,:,2)=rduc(:,:,:,2)+coeff(:,:,:,21,step)*rvtemp(:,:,:,1) &
00347           &                   +coeff(:,:,:,22,step)*rvtemp(:,:,:,2) &
00348           &                   +coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00349      ! - (dr . nabla ) U_0
00350      rduc(:,:,:,2)=rduc(:,:,:,2)-rd(:,:,:,1)*coeff(:,:,:,25,nonstep)&
00351           &                   -rd(:,:,:,2)*coeff(:,:,:,28,nonstep)&
00352           &                   -rd(:,:,:,3)*coeff(:,:,:,31,nonstep)
00353      !end y_component
00354      !third z_component
00355      ! + (U_0 . nabla ) dr
00356      !calculate d/dx_i rd_j  and move into real space
00357      fvtemp(:,:,:,1)=ddxe(fd(:,:,:,3))
00358      fvtemp(:,:,:,2)=ddye(fd(:,:,:,3))
00359      call sfftw_execute(plan_f2rvtemp_1)
00360      call sfftw_execute(plan_f2rvtemp_2)
00361      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00362      rvtemp(:,:,:,3)=ddze(rd(:,:,:,3))
00363      rduc(:,:,:,3)=rduc(:,:,:,3)+coeff(:,:,:,21,step)*rvtemp(:,:,:,1) &
00364           &                     +coeff(:,:,:,22,step)*rvtemp(:,:,:,2) &
00365           &                     +coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00366      ! - (dr . nabla ) U_0
00367      rduc(:,:,:,3)=rduc(:,:,:,3)-rd(:,:,:,1)*coeff(:,:,:,26,nonstep)&
00368           &                     -rd(:,:,:,2)*coeff(:,:,:,29,nonstep)&
00369           &                     -rd(:,:,:,3)*coeff(:,:,:,32,nonstep)
00370      !end z_component
00371 
00372 
00373      !======================================================================
00374      ! now cell ends 
00375      rdUe(:,:,:,:)=rd(:,:,:,4:6)
00376      !first x_component
00377      ! + (U_0 . nabla ) dr
00378      !calculate d/dx_i rd_j  and move into real space
00379      fvtemp(:,:,:,1)=ddxe(fd(:,:,:,1))
00380      fvtemp(:,:,:,2)=ddye(fd(:,:,:,1))
00381      call sfftw_execute(plan_f2rvtemp_1)
00382      call sfftw_execute(plan_f2rvtemp_2)
00383      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00384      rvtemp(:,:,:,3)=ddze(rd(:,:,:,1))
00385      rdue(:,:,:,1)=rdue(:,:,:,1)+coeff(:,:,:,21,nonstep)*rvtemp(:,:,:,1) &
00386           &                     +coeff(:,:,:,22,nonstep)*rvtemp(:,:,:,2) &
00387           &                     +coeff(:,:,:,23,nonstep)*rvtemp(:,:,:,3) 
00388      ! - (dr . nabla ) U_0
00389      rdue(:,:,:,1)=rdue(:,:,:,1)-rd(:,:,:,1)*coeff(:,:,:,24,nonstep)&
00390           &                     -rd(:,:,:,2)*coeff(:,:,:,27,nonstep)&
00391           &                     -rd(:,:,:,3)*coeff(:,:,:,30,nonstep)
00392      !end x_component
00393      ! second y_component
00394      ! + (U_0 . nabla ) dr
00395      !calculate d/dx_i rd_j  and move into real space
00396      fvtemp(:,:,:,1)=ddxe(fd(:,:,:,2))
00397      fvtemp(:,:,:,2)=ddye(fd(:,:,:,2))
00398      call sfftw_execute(plan_f2rvtemp_1)
00399      call sfftw_execute(plan_f2rvtemp_2)
00400      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00401      rvtemp(:,:,:,3)=ddze(rd(:,:,:,2))
00402      rdue(:,:,:,2)=rdue(:,:,:,2)+coeff(:,:,:,21,nonstep)*rvtemp(:,:,:,1) &
00403           &                     +coeff(:,:,:,22,nonstep)*rvtemp(:,:,:,2) &
00404           &                     +coeff(:,:,:,23,nonstep)*rvtemp(:,:,:,3) 
00405      ! - (dr . nabla ) U_0
00406      rdue(:,:,:,2)=rdue(:,:,:,2)-rd(:,:,:,1)*coeff(:,:,:,25,nonstep)&
00407           &                     -rd(:,:,:,2)*coeff(:,:,:,28,nonstep)&
00408           &                     -rd(:,:,:,3)*coeff(:,:,:,31,nonstep)
00409      !end y_component
00410      !third z_component
00411      ! + (U_0 . nabla ) dr
00412      !calculate d/dx_i rd_j  and move into real space
00413      fvtemp(:,:,:,1)=ddxe(fd(:,:,:,3))
00414      fvtemp(:,:,:,2)=ddye(fd(:,:,:,3))
00415      call sfftw_execute(plan_f2rvtemp_1)
00416      call sfftw_execute(plan_f2rvtemp_2)
00417      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00418      rvtemp(:,:,:,3)=ddze(rd(:,:,:,3))
00419      rdue(:,:,:,3)=rdue(:,:,:,3)+coeff(:,:,:,21,nonstep)*rvtemp(:,:,:,1) &
00420           &                     +coeff(:,:,:,22,nonstep)*rvtemp(:,:,:,2) &
00421           &                     +coeff(:,:,:,23,nonstep)*rvtemp(:,:,:,3) 
00422      ! - (dr . nabla ) U_0
00423      rdue(:,:,:,3)=rdue(:,:,:,3)-rd(:,:,:,1)*coeff(:,:,:,26,nonstep)&
00424           &                     -rd(:,:,:,2)*coeff(:,:,:,29,nonstep)&
00425           &                     -rd(:,:,:,3)*coeff(:,:,:,32,nonstep)
00426      !end z_component
00427 
00428      !=================We now have dU (rs, centred and cell ends)=======
00429      !BC on dU
00430      !(U_0 assumed to vanish at vertical boundaries) => bc's OK
00431 
00432 
00433      !===========================BCs=======================================
00434      ! SHOULD BE NO NEED - WE ASSUME U_0=0 at Boundaries
00435 
00436      !copy into fourier space as well 
00437      call sfftw_execute(plan_r2fdUc_1)
00438      call sfftw_execute(plan_r2fdUc_2)
00439      call sfftw_execute(plan_r2fdUc_3)
00440      call sfftw_execute(plan_r2fdUe_1)
00441      call sfftw_execute(plan_r2fdUe_2)
00442      call sfftw_execute(plan_r2fdUe_3)
00443 
00444      !---------------------------------------------------------------------------
00445      !----------------------Calculate the force due to U_0-----------------------
00446 
00447      !rinertia=[ -(U_0 . nabla) d/dt(dr) - (U_0 . nabla) dU 
00448      !           +([d/dt(dr) -drho U_0/rho_0 -dU].nabla) U_0]
00449 
00450      !x-component
00451      !-(U_0 . nabla) (d/dt(dr))
00452      fvtemp(:,:,:,1)=ddxc(fd(:,:,:,4))
00453      fvtemp(:,:,:,2)=ddyc(fd(:,:,:,4))
00454      call sfftw_execute(plan_f2rvtemp_1)
00455      call sfftw_execute(plan_f2rvtemp_2)
00456      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00457      rvtemp(:,:,:,3)=ddzc(rd(:,:,:,4))
00458      rinertia(:,:,:,1)=-coeff(:,:,:,21,step)*rvtemp(:,:,:,1)&
00459           &            -coeff(:,:,:,22,step)*rvtemp(:,:,:,2) &
00460           &            -coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00461      !-(U_0 . nabla) dU
00462      fvtemp(:,:,:,1)=ddxc(fdUe(:,:,:,1))
00463      fvtemp(:,:,:,2)=ddyc(fdUe(:,:,:,1))
00464      call sfftw_execute(plan_f2rvtemp_1)
00465      call sfftw_execute(plan_f2rvtemp_2)
00466      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00467      rvtemp(:,:,:,3)=ddzc(rdUe(:,:,:,1))
00468      rinertia(:,:,:,1)=rinertia(:,:,:,1)-coeff(:,:,:,21,step)*rvtemp(:,:,:,1)&
00469           &                             -coeff(:,:,:,22,step)*rvtemp(:,:,:,2)&
00470           &                             -coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00471 
00472      !+(d/dt(dr) . nabla) U_0]
00473 
00474      rinertia(:,:,:,1)=rinertia(:,:,:,1)&
00475           &               +centre(rd(:,:,:,4))*coeff(:,:,:,24,nonstep)&
00476           &               +centre(rd(:,:,:,5))*coeff(:,:,:,27,nonstep)&
00477           &               +centre(rd(:,:,:,6))*coeff(:,:,:,30,nonstep)
00478 
00479      !+([ -drho U_0/rho_0 -dU].nabla) U_0]
00480      rvtemp(:,:,:,1)=-rrhoc*coeff(:,:,:,33,step)-rdUc(:,:,:,1)
00481      rvtemp(:,:,:,2)=-rrhoc*coeff(:,:,:,34,step)-rdUc(:,:,:,2)
00482      rvtemp(:,:,:,3)=-rrhoc*coeff(:,:,:,35,step)-rdUc(:,:,:,3)
00483 
00484      rinertia(:,:,:,1)=rinertia(:,:,:,1)+rvtemp(:,:,:,1)*coeff(:,:,:,24,step)&
00485           &                             +rvtemp(:,:,:,2)*coeff(:,:,:,27,step)&
00486           &                             +rvtemp(:,:,:,3)*coeff(:,:,:,30,step)
00487      !------------------------END X-----------------------------------------
00488      !y-component
00489      !-(U_0 . nabla) d/dt(dr)
00490      fvtemp(:,:,:,1)=ddxc(fd(:,:,:,5))
00491      fvtemp(:,:,:,2)=ddyc(fd(:,:,:,5))
00492      call sfftw_execute(plan_f2rvtemp_1)
00493      call sfftw_execute(plan_f2rvtemp_2)
00494      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00495      rvtemp(:,:,:,3)=ddzc(rd(:,:,:,5))
00496      rinertia(:,:,:,2)=-coeff(:,:,:,21,step)*rvtemp(:,:,:,1)&
00497           &            -coeff(:,:,:,22,step)*rvtemp(:,:,:,2) &
00498           &            -coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00499      !-(U_0 . nabla) dU
00500      fvtemp(:,:,:,1)=ddxc(fdUe(:,:,:,2))
00501      fvtemp(:,:,:,2)=ddyc(fdUe(:,:,:,2))
00502      call sfftw_execute(plan_f2rvtemp_1)
00503      call sfftw_execute(plan_f2rvtemp_2)
00504      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00505      rvtemp(:,:,:,3)=ddzc(rdUe(:,:,:,2))
00506      rinertia(:,:,:,2)=rinertia(:,:,:,2)-coeff(:,:,:,21,step)*rvtemp(:,:,:,1)&
00507           &                             -coeff(:,:,:,22,step)*rvtemp(:,:,:,2)&
00508           &                             -coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00509      !+(d/dt(dr) . nabla) U_0]
00510      rinertia(:,:,:,2)=rinertia(:,:,:,2)&
00511           &                   +centre(rd(:,:,:,4))*coeff(:,:,:,25,nonstep)&
00512           &                   +centre(rd(:,:,:,5))*coeff(:,:,:,28,nonstep)&
00513           &                   +centre(rd(:,:,:,6))*coeff(:,:,:,31,nonstep)
00514 
00515      !+([-drho U_0/rho_0 -dU].nabla) U_0]
00516      rvtemp(:,:,:,1)=-rrhoc*coeff(:,:,:,33,step)-rdUc(:,:,:,1)
00517      rvtemp(:,:,:,2)=-rrhoc*coeff(:,:,:,34,step)-rdUc(:,:,:,2)
00518      rvtemp(:,:,:,3)=-rrhoc*coeff(:,:,:,35,step)-rdUc(:,:,:,3)
00519      rinertia(:,:,:,2)=rinertia(:,:,:,2)+rvtemp(:,:,:,1)*coeff(:,:,:,25,step)&
00520           &                             +rvtemp(:,:,:,2)*coeff(:,:,:,28,step)&
00521           &                             +rvtemp(:,:,:,3)*coeff(:,:,:,31,step)
00522      !------------------------END y-----------------------------------------
00523      !z-component
00524      !-(U_0 . nabla) d/dt(dr)
00525      fvtemp(:,:,:,1)=ddxc(fd(:,:,:,6))
00526      fvtemp(:,:,:,2)=ddyc(fd(:,:,:,6))
00527      call sfftw_execute(plan_f2rvtemp_1)
00528      call sfftw_execute(plan_f2rvtemp_2)
00529      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00530      rvtemp(:,:,:,3)=ddzc(rd(:,:,:,6))
00531 
00532      rinertia(:,:,:,3)=-coeff(:,:,:,21,step)*rvtemp(:,:,:,1)&
00533           &            -coeff(:,:,:,22,step)*rvtemp(:,:,:,2)&
00534           &            -coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00535      !-(U_0 . nabla) dU
00536      fvtemp(:,:,:,1)=ddxc(fdUe(:,:,:,3))
00537      fvtemp(:,:,:,2)=ddyc(fdUe(:,:,:,3))
00538      call sfftw_execute(plan_f2rvtemp_1)
00539      call sfftw_execute(plan_f2rvtemp_2)
00540      rvtemp(:,:,:,1:2)=rvtemp(:,:,:,1:2)*norm
00541      rvtemp(:,:,:,3)=ddzc(rdUe(:,:,:,3))
00542      rinertia(:,:,:,3)=rinertia(:,:,:,3)-coeff(:,:,:,21,step)*rvtemp(:,:,:,1)&
00543           &                             -coeff(:,:,:,22,step)*rvtemp(:,:,:,2)&
00544           &                             -coeff(:,:,:,23,step)*rvtemp(:,:,:,3) 
00545      !+(d/dt(dr) . nabla) U_0]
00546      rinertia(:,:,:,3)=rinertia(:,:,:,3)&
00547           &                     +centre(rd(:,:,:,4))*coeff(:,:,:,26,nonstep)&
00548           &                     +centre(rd(:,:,:,5))*coeff(:,:,:,29,nonstep)&
00549           &                     +centre(rd(:,:,:,6))*coeff(:,:,:,32,nonstep)    
00550 
00551      !+([-drho U_0/rho_0 -dU].nabla) U_0]
00552      rvtemp(:,:,:,1)=-rrhoc*coeff(:,:,:,33,step)-rdUc(:,:,:,1)
00553      rvtemp(:,:,:,2)=-rrhoc*coeff(:,:,:,34,step)-rdUc(:,:,:,2)
00554      rvtemp(:,:,:,3)=-rrhoc*coeff(:,:,:,35,step)-rdUc(:,:,:,3)
00555      rinertia(:,:,:,3)=rinertia(:,:,:,3)+rvtemp(:,:,:,1)*coeff(:,:,:,26,step)&
00556           &                             +rvtemp(:,:,:,2)*coeff(:,:,:,29,step)&
00557           &                             +rvtemp(:,:,:,3)*coeff(:,:,:,32,step)  
00558      !------------------------END z-----------------------------------------
00559      !BC on forces?
00560 
00561 
00562      !========================END U_0 force stuff==============================
00563 
00564   endif !!!!!!!!!!END FLOW!!!!!!!!!!!!!!!!!!!!
00565 
00566 
00567 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00568 !!!!!!!!!!!NEXT FIND d/dt (drdot)
00569   !-d/dx (P) :FS
00570   fvtemp(:,:,:,1)=-ddxc(fpress)
00571   !-d/dy (P) :FS
00572   fvtemp(:,:,:,2)=-ddyc(fPress)       
00573   !-d/dz (P) :RS
00574   rvtemp(:,:,:,3)=-ddzc(rpress)
00575   !d/dx (P), d/dy (P) ->RS
00576   call sfftw_execute(plan_f2rvtemp_1)
00577   call sfftw_execute(plan_f2rvtemp_2)
00578   rvtemp(:,:,:,1)=rvtemp(:,:,:,1)*norm
00579   rvtemp(:,:,:,2)=rvtemp(:,:,:,2)*norm
00580   drddt(:,:,:,4)=rvtemp(:,:,:,1)*coeff(:,:,:,12,step)  &
00581        &     +rLor(:,:,:,1)*coeff(:,:,:,12,nonstep)            &
00582        &     +rInertia(:,:,:,1)
00583   drddt(:,:,:,5)=rvtemp(:,:,:,2)*coeff(:,:,:,12,step)  &
00584        &     +rLor(:,:,:,2)*coeff(:,:,:,12,nonstep)            &
00585        &     +rinertia(:,:,:,2)
00586   drddt(:,:,:,6)=rvtemp(:,:,:,3)*coeff(:,:,:,12,step)  &
00587        &     +rLor(:,:,:,3)*coeff(:,:,:,12,step)            &
00588        &     +rinertia(:,:,:,3)                        &
00589        &          +rrhoc*coeff(:,:,:,13,step)
00590 
00591   !ENDBC
00592 
00593   call sfftw_execute(plan_r2fddt_4)
00594   call sfftw_execute(plan_r2fddt_5)
00595   call sfftw_execute(plan_r2fddt_6)
00596 
00597 
00598 end subroutine derivs