molecular and hyperviscosity

src/makefile

@@ -35,7 +35,6 @@ timestep: $(objects) $(sourcedir)/timestep.f90
 bminmax: $(objects) $(sourcedir)/bminmax.f90
 	$(f90) parameters.o constants.o $(sourcedir)/bminmax.f90 -o bminmax $(flags)
 
-
 $(objects): $(fft_lib) $(sources)
 	$(f90) $(fft_lib) $(sources) -c $(flags)
 

src/parameters.f90

@@ -4,20 +4,27 @@ module parameters
 ! the suite of ps f90 files.
 
  !Domain grid dimensions:
-integer,parameter:: nx=256,ny=384
+integer,parameter:: nx=N_X,ny=N_Y
 
  !Domain width in x and limits in y:
-double precision,parameter:: ellx=1000d0
-double precision,parameter:: ymin=0.d0,ymax=1500d0
+double precision,parameter:: ellx=L_X
+double precision,parameter:: ymin=Y_MIN,ymax=Y_MAX
 
  !Simulation duration and data save interval:
-double precision,parameter:: tsim=900.d0,tgsave=10.d0
+double precision,parameter:: tsim=T_SIM,tgsave=T_GSAVE
 
- !Hyperviscosity parameters:
-integer,parameter:: nnu=3
-double precision,parameter:: prediss=10.d0
-! Damping rate is prediss*zeta_rms*(k/k_max)^(2*nnu) on wavenumber
-! k where k_max is the maximum x or y wavenumber and zeta_rms is
-! the rms vorticity. Note: nnu = 3 and prediss = 10 are recommended.
+ !(Hyper)viscosity parameters:
+integer,parameter:: nnu=N_NU
+double precision,parameter:: prediss=PRE_DISS
+! If nnu = 1, this is the molecular viscosity case.  Then, we 
+! choose the viscosity nu = prediss*((b_max-b_min)/k_{x,max}^3)
+! where k_{x_max} is the maximum x wavenumber.
+! Note: prediss = 2 is recommended.
+! ----------------------------------------------------------------
+! If nnu > 1, this is the hyperviscosity case.  Then, the damping
+! rate is prediss*zeta_char*(k/k_max)^(2*nnu) on wavenumber k
+! where k_max is the maximum x or y wavenumber and zeta_char is
+! a characteristic vorticity (see subroutine adapt of strat.f90).
+! Note: nnu = 3 and prediss = 10 are recommended.
 
 end module

src/spectral.f90

@@ -17,19 +17,22 @@ module spectral
 integer:: kmag(0:nxm1,0:ny),kmax
 
 !====================================================================!
-! From main code: call init_spectral          to initialise          !
+! From main code: call init_spectral(bbdif)   to initialise          !
 ! then            call main_invert(zz,uu,vv)  to perform inversion   !
 !====================================================================!
 
 contains
 
 !=====================================================================
 
-subroutine init_spectral
+subroutine init_spectral(bbdif)
 
  !Declarations:
 implicit none
 
+ !Passed variable (bbdif = max(b) - min(b) at t = 0):
+double precision:: bbdif
+
  !Local variables:
 
  !Maximum x wavenumber:
@@ -154,20 +157,45 @@ subroutine init_spectral
 enddo
 
 !---------------------------------------------------------------------
- !Define dissipation operator:
-visc=prediss/max(rkxmax,rkymax)**(2*nnu)
-hdis(0,0)=zero
-do kx=1,nxm1
-  hdis(kx,0)=visc*rkx(kx)**(2*nnu)
-enddo
-do ky=1,ny
-  hdis(0,ky)=visc*rky(ky)**(2*nnu)
-enddo
-do ky=1,ny
+! Damping, viscous or hyperviscous:
+if (nnu .eq. 1) then
+   !Define viscosity:
+  visc=prediss*sqrt(bbdif/rkxmax**3)
+  write(*,'(a,1p,e14.7)') ' Viscosity nu = ',visc
+
+   !Define spectral dissipation operator:
+  hdis(0,0)=zero
   do kx=1,nxm1
-    hdis(kx,ky)=visc*(rkx(kx)**2+rky(ky)**2)**nnu
+    hdis(kx,0)=visc*rkx(kx)**2
   enddo
-enddo
+  do ky=1,ny
+    hdis(0,ky)=visc*rky(ky)**2
+  enddo
+  do ky=1,ny
+    do kx=1,nxm1
+      hdis(kx,ky)=visc*(rkx(kx)**2+rky(ky)**2)
+    enddo
+  enddo
+
+else
+   !Define hyperviscosity:
+  visc=prediss/max(rkxmax,rkymax)**(2*nnu)
+  write(*,'(a,1p,e14.7)') ' Hyperviscosity nu = ',visc
+
+   !Define dissipation operator:
+  hdis(0,0)=zero
+  do kx=1,nxm1
+    hdis(kx,0)=visc*rkx(kx)**(2*nnu)
+  enddo
+  do ky=1,ny
+    hdis(0,ky)=visc*rky(ky)**(2*nnu)
+  enddo
+  do ky=1,ny
+    do kx=1,nxm1
+      hdis(kx,ky)=visc*(rkx(kx)**2+rky(ky)**2)**nnu
+    enddo
+  enddo
+endif
 
 return 
 end subroutine init_spectral

src/strat.f90

@@ -93,31 +93,33 @@ subroutine initialise
 
 implicit none
 
- !Local variable:
+ !Local variables:
 double precision:: qq(0:ny,0:nxm1)
+double precision:: bbdif
 
 !----------------------------------------------------------------------
- !Initialise inversion constants and arrays:
-call init_spectral
-
-!----------------------------------------------------------------------
- !Read vorticity and convert to spectral space as zs:
-open(11,file='zz_init.r8',form='unformatted', &
+ !Read buoyancy and find min/max value for use in determining viscosity:
+open(11,file='bb_init.r8',form='unformatted', &
     & access='direct',status='old',recl=2*nbytes)
 read(11,rec=1) t,qq
 close(11)
 
-call ptospc_fc(nx,ny,qq,zs,xfactors,yfactors,xtrig,ytrig)
+bbdif=maxval(qq)-minval(qq)
+
+ !Initialise inversion constants and arrays:
+call init_spectral(bbdif)
+
+ !Convert buoyancy (in qq) to spectral space as bs:
+call ptospc_fc(nx,ny,qq,bs,xfactors,yfactors,xtrig,ytrig)
  !qq is overwritten here.
 
-!----------------------------------------------------------------------
- !Read buoyancy and convert to spectral space as bs:
-open(11,file='bb_init.r8',form='unformatted', &
+ !Read vorticity and convert to spectral space as zs:
+open(11,file='zz_init.r8',form='unformatted', &
     & access='direct',status='old',recl=2*nbytes)
 read(11,rec=1) t,qq
 close(11)
 
-call ptospc_fc(nx,ny,qq,bs,xfactors,yfactors,xtrig,ytrig)
+call ptospc_fc(nx,ny,qq,zs,xfactors,yfactors,xtrig,ytrig)
  !qq is overwritten here.
 
 !----------------------------------------------------------------------
@@ -383,13 +385,21 @@ subroutine adapt
 dt4=dt/four
 
 !---------------------------------------------------------------------
- !Update hyperdiffusion operator used in time stepping:
-dfac=zzch*dt/two
- !zzch is the characteristic vorticity defined above.
-diss=two/(one+dfac*hdis)
- !hdis = C*(K/K_max)^{2p} where K^2 = k_x^2+k_y^2, p is the order,
- !K_max is the maximum x or y wavenumber and C is a dimensionless
- !prefactor (see spectral.f90 and parameters.f90 where C = prediss).
+if (nnu .eq. 1) then
+  !Update diffusion operator used in time stepping:
+  dfac=dt/two
+  diss=two/(one+dfac*hdis)
+   !hdis = nu*(k_x^2+k_y^2) where nu is the viscosity coefficient
+   !(see spectral.f90 and parameters.f90).
+else
+   !Update hyperdiffusion operator used in time stepping:
+  dfac=zzch*dt/two
+   !zzch is the characteristic vorticity defined above.
+  diss=two/(one+dfac*hdis)
+   !hdis = C*(K/K_max)^{2p} where K^2 = k_x^2+k_y^2, p is the order,
+   !K_max is the maximum x or y wavenumber and C is a dimensionless
+   !prefactor (see spectral.f90 and parameters.f90 where C = prediss).
+endif
 
 !---------------------------------------------------------------------
  !Save |u|_max, N_max and gamma_max to monitor.asc: