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

00001 !================================================================
00002 subroutine derivs(PE,KE)
00003   ! Calculates RHS
00004   implicit none     
00005   include 'param.f90'
00006   include 'terms3d.f90'           
00007   complex fdtemp(nx/2+1,ny,nz,nvar)
00008   real  PE,KE
00009   integer i,k,n
00010 !FFT IS NOT NORMALIZED  
00011 !!!!FIRST FIND THE EULERIAN DENSITY PERTURBATION!!!!!!!!!!!!!!!
00012   fdtemp=fd
00013   ! dr -> RS
00014 !bottom bc is no displacement
00015   call sfftw_execute(plan_f2rd_1)
00016   call sfftw_execute(plan_f2rd_2)
00017   call sfftw_execute(plan_f2rd_3)
00018   call sfftw_execute(plan_f2rd_4)
00019   call sfftw_execute(plan_f2rd_5)
00020   call sfftw_execute(plan_f2rd_6)
00021   fd=fdtemp
00022   rd=rd*norm
00023 ! move to variables not scaled by rho_0
00024   do i=1,6 
00025      rd(:,:,:,i)=rd(:,:,:,i)*coeff(:,:,:,12)
00026   enddo
00027  
00028   PE=sum(abs(rd(:,:,:,1)))
00029   KE=sum(abs(rd(:,:,:,4)))
00030 
00031   ! rho_0 dr :RS
00032   rvtemp(:,:,:,1)=coeff(:,:,:,1)*rd(:,:,:,1)
00033   rvtemp(:,:,:,2)=coeff(:,:,:,1)*rd(:,:,:,2)
00034   rvtemp(:,:,:,3)=coeff(:,:,:,1)*rd(:,:,:,3)
00035   ! d/dz(rho_0 dr_z) :RS
00036   rvtemp(:,:,:,3)=ddz(rvtemp(:,:,:,3))
00037   ! rho_0 dr_x, rho_0 dr_y -> FS
00038   call sfftw_execute(plan_r2fvtemp_1)
00039   call sfftw_execute(plan_r2fvtemp_2)
00040 !  call sfftw_execute(plan_r2fvtemp_3)
00041 !!!!! NOW STORE FD/dens in fdtemp 
00042 !  fdtemp=fvtemp 
00043   ! d/dx(rho_0 dr_x) :FS
00044   fvtemp(:,:,:,1)=ddx(fvtemp(:,:,:,1))
00045   ! d/dy(rho_0 dr_y) :FS       
00046   fvtemp(:,:,:,2)=ddy(fvtemp(:,:,:,2))
00047   ! d/dx(rho_0 dr_x) ->RS
00048   call sfftw_execute(plan_f2rvtemp_1)
00049   rvtemp(:,:,:,1)=rvtemp(:,:,:,1)*norm
00050   ! d/dz(rho_0 dr_z) ->RS
00051   call sfftw_execute(plan_f2rvtemp_2)
00052   rvtemp(:,:,:,2)=rvtemp(:,:,:,2)*norm
00053   rrho=-(rvtemp(:,:,:,1)+rvtemp(:,:,:,2)+rvtemp(:,:,:,3))
00054 !------------------BC-------------------: RS
00055   rrho(:,:,nz)=-( rd(:,:,nz,1)*coeff(:,:,nz,2) &
00056                  +rd(:,:,nz,2)*coeff(:,:,nz,3) &
00057                  +rd(:,:,nz,3)*coeff(:,:,nz,4))
00058  call sfftw_execute(plan_r2frho)
00059         do n=1,nx/2+1
00060            frho(n,:,1)=frho(n,:,2)*bc(n)
00061         enddo
00062  call sfftw_execute(plan_f2rrho)
00063   rrho(:,:,1)=rrho(:,:,1)*coeff(:,:,1,4)/coeff(:,:,2,4)
00064   rrho=rrho*norm
00065 ! ENDBC
00066 
00067 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00068 !!!!!!NEXT FIND THE EULERIAN PRESSURE PERTURBATION!!!!!!!!!!!!
00069 !(GAMMA_1 P_0/rho_0) rho_1 + (GAMMA_1 P_0 Grad(rho_0)/rho_0 - Grad(P_0)) dr :RS
00070   rPress=coeff(:,:,:,5)*rrho &
00071         +coeff(:,:,:,6)*rd(:,:,:,1)&
00072         +coeff(:,:,:,7)*rd(:,:,:,2)&
00073         +coeff(:,:,:,8)*rd(:,:,:,3)
00074 !------------------BC-------------------:RS
00075 !
00076   rpress(:,:,nz)=-( rd(:,:,nz,1)*coeff(:,:,nz,9) &
00077                    +rd(:,:,nz,2)*coeff(:,:,nz,10) &
00078                    +rd(:,:,nz,3)*coeff(:,:,nz,11)) 
00079  call sfftw_execute(plan_r2fpress)
00080         do n=1,nx/2+1
00081            fpress(n,:,1)=fpress(n,:,2)*bc(n)
00082         enddo  
00083 call sfftw_execute(plan_f2rpress)
00084   rpress(:,:,1)=rpress(:,:,1)*coeff(:,:,1,11)/coeff(:,:,2,11)
00085   rpress=rpress*norm
00086 ! ENDBC
00087 
00088 !rpress ->FS       
00089 !   call sfftw_execute(plan_r2fpress)
00090 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00091 !!!!!!!!!!!NEXT FIND d/dt (drdot)
00092 !d/dx (P) :FS
00093    fvtemp(:,:,:,1)=ddx(fpress)
00094 !d/dy (P) :FS
00095    fvtemp(:,:,:,2)=ddy(fPress)       
00096 !d/dz (P) :RS
00097    rvtemp(:,:,:,3)=ddz(rpress)
00098 !d/dx (P), d/dy (P) ->RS
00099   call sfftw_execute(plan_f2rvtemp_1)
00100   call sfftw_execute(plan_f2rvtemp_2)
00101   rvtemp(:,:,:,1)=rvtemp(:,:,:,1)*norm
00102   rvtemp(:,:,:,2)=rvtemp(:,:,:,2)*norm
00103 ! Grad(P)/rho_0 +  rho g/rho_0 (hat(y))
00104 
00105   drddt(:,:,:,1)=rd(:,:,:,4)
00106   drddt(:,:,:,2)=rd(:,:,:,5)
00107   drddt(:,:,:,3)=rd(:,:,:,6)
00108   drddt(:,:,:,4)=-rvtemp(:,:,:,1)*coeff(:,:,:,12)
00109   drddt(:,:,:,5)=-rvtemp(:,:,:,2)*coeff(:,:,:,12)
00110   drddt(:,:,:,6)=-rvtemp(:,:,:,3)*coeff(:,:,:,12) + rrho*coeff(:,:,:,13)
00111 
00112 !move to variables scaled by rho_0
00113   do i=1,6
00114     drddt(:,:,:,i)=drddt(:,:,:,i)*coeff(:,:,:,1)
00115   enddo 
00116   call sfftw_execute(plan_r2fddt_1)
00117   call sfftw_execute(plan_r2fddt_2)
00118   call sfftw_execute(plan_r2fddt_3)   
00119   call sfftw_execute(plan_r2fddt_4)
00120   call sfftw_execute(plan_r2fddt_5)
00121   call sfftw_execute(plan_r2fddt_6)
00122 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00123 !!!!!!!LASTLY d/dt (dr)
00124 !  dfddt(:,:,:,1:3)=fd(:,:,:,4:6)
00125 
00126 end subroutine derivs