source: trunk/SpectralModelF90/PROJECTS/SIMU_ASARO/TEST2/user_module.f90 @ 2

module user_module    
!! Routines to prescribe ics and forcing 
!! as functions of space/time
!! flow_solve will include user_module.
 public  :: user_defined_ics         
 public  :: user_defined_forcing 
    
!!defns of background scalar profiles
!!s1_bar(z), s2_bar(z)
 public  :: eval_s1_bar
 public  :: eval_s2_bar
        
 contains

!===========================================================================
    subroutine user_defined_ics(x,y,z,Lx,Ly,Lz, &
                                u,v,w,s1,s2,f,beta,U0,  &
                                g,rho_0,DGRAD)
!
!     inputs 
!      x,y,z     position in spatial domain            [m]
!      Lx,Ly,Lz  domain lengths in x,y,z               [m]
!      f         Coriolis param (set in problem_params)[1/s]
!      g         gravity                 ""            [m/s2]
!      rho_0     reference fluid density ""            [kg/m3]
!      DGRAD     characteristic density gradient  ""   [kg/m4]
!
!     outputs
!      u,v,w     speeds in x,y,z directions at (x,y,z)     [m/s] 
!      s1        e.g. could be perturbation temperature in [deg C]
!      s2        e.g. could be perturbation salinity in    [psu]
!
!     note:
!      s1 and s2 are defined relative to the ambient profiles
!      s1_bar and s2_bar defined in the routines eval_s1_bar
!      and eval_s2_bar below. There are corresponding quantities
!      rho_bar and perturbation density pd that are computed
!      via the eqn of state (routine eos) and used internally
!      to compute buoyancy forces. They need not be set here.

      implicit none
      real(kind=8)          :: x,y,z,Lx,Ly,Lz,u,v,w,s1,s2
      real(kind=8)          :: f,beta,g,rho_0,DGRAD,pi,fx,U0
!===================================================================
!       user declares any additional LOCAL variables needed here
      real(kind=8)                      :: arg,c11,c12,c21,c22,omega_ics
      real(kind=8)                      :: A_ics,k_ics,l_ics,m_ics,N2_ics,omega2

!===================================================================
! Example:
! Initial conditions are a snapshot of an internal wave
! in a uniform density stratification. I've got the scalars
! set to "density" and "passive tracer". The ambient density
! gradient (actually its absolute value) is a user specified
! parameter in problem_params. Its value is made available
! here in case it is needed. (Also f,g,rho_0) 

       u=0.
       v=0.
       w=0.
       s1=0.
       s2=0.

       return
      end subroutine user_defined_ics

      
 !===============================================================================     
      subroutine user_defined_forcing(x,y,z,t,F1,F2,F3,F4,F5,  &
                                      Lx,Ly,Lz,f,istep,iend_storm)
!     User defined subroutine to add time/space dependent forcing to the 
!     equations of motion. Inputs and outputs for this routine are all in
!     dimensional units.
!

!     inputs:
!      x,y,z     spatial location at which forcing is to be defined, [m]
!      t         time at which forcing is to be defined              [s]
!      f,g,rho_0 Coriolis parameter, gravity, reference density [1/s,m/s2,kg/m3]

!     outputs:
!      F1-3      acceleration at x,y,z & t in x,y,z directions       [m/s2]
!      F4        rhs of perturbation scalar (s1) eqn                 [(deg C)/s]
!      F5        rhs of perturbation scalar (s2) eqn                 [(psu)/s]
!                careful with F4,F5 you must take into account what you are
!                using for eval_s1_bar & eval_s2_bar when you define this term
!===============================================================================

      implicit none
      integer           ::  istep,iend_storm
      real(kind=8)      ::  x,y,z
      real(kind=8)      ::  u,v,w,s1,s2,F1,F2,F3,F4,F5
      real(kind=8)      ::  t,Lx,Ly,Lz,f,g,rho_0

!      real(kind=8),dimension(:),intent(in)      ::  x,y,z
!      real(kind=8),dimension(:,:),intent(in)      ::  u,v,w,s1,s2
!      real(kind=8)      ::  Lx,Ly,Lz,f,g,rho_0

!!      local variable
      real(kind=8)      ::   A_ics,xfreq,arg,sigma1,sigma2,FF,GG
 
      if ( istep .le. iend_storm ) then

!=================FORCAGE XLV=======  
!#      if ( x == 0. .and. z == 0. .and. y == 0. ) then
        print*,'========='
        print*,' Forcage, istep=',istep
        print*,'========='
!#      endif

       A_ics=100.                                     ! displacement [m]
       xfreq=1.2
       arg= xfreq*f*t  ! generally -omega_ics*t but set t=0,
       sigma1=20000.
       sigma2=100.
       FF=exp((-((x+Lx/2)**2))/(sigma1**2))
       GG=exp((-((Lz/2-z)**2))/(sigma2**2))

!     Define rhs values for evolution equations
       s1=0.
       s2=0.
       w = 0.
       u = -A_ics*xfreq*xfreq*f*f*sin(arg)*FF*GG          ! du/dt
       v = 0.
!===================================

      else

      F1= 0.
      F2= 0.
      F3= 0.
      F4= 0.
      F5= 0.

      endif

      return
     end subroutine user_defined_forcing
!============================================================================


!============================================================================
subroutine initialize_particles( positions,npts,Lx,Ly,Lz )
     use Ziggurat
    
     implicit none
     integer,  parameter  ::  DP=SELECTED_REAL_KIND( 12, 60 )
     integer        :: npts,i,j            !! the number of particles
     integer        :: myid,ierr 
     real(kind=8)   :: positions(npts,3)   !! initial particle positions [m]
     real(kind=8)   :: Lx,Ly,Lz            !! domain size in [m]
     real(DP)       :: urv(3),L(3)         !! uniform random deviates in [0,1]
     include 'mpif.h'
     
     !!   call zigset( iseed ) executed once during initialization
     !!   rnor( ),rexp() also available
     
     call mpi_comm_rank (MPI_COMM_WORLD,myid,ierr)
     
     L(:)=(/Lx,Ly,Lz/)
     do i=1,npts
      do j=1,3
       urv(j) = uni( )
       positions(i,j)=L(j)*urv(j)
      enddo
     if(myid==0) write(0,*) i,positions(i,:)/Lz
     enddo
     
   end subroutine initialize_particles
!============================================================================   
   
   
   
!=============================================================================     
      subroutine eval_s1_bar(z,Lz,s1_bar,s1_bar_z,s1_bar_zz)  
!      input
!      z            vertical z spatial position  [m]
!      Lz           vertical size of computational domain [m]
!
!     output 
!      s1_bar, s1_bar_z, s1_bar_zz
!      scalar concentration and z derivs   [dimensional units]

      implicit none
      real(kind=8) s1_bar,s1_bar_z,s1_bar_zz,z,Lz,dTdz,s1_min

!     LOCAL variables: you may need some intermediate variables,
!                      these need to be explicitly declared here
      real(kind=8)   :: DTdz = 2.1e-3
      real(kind=8)   :: s1_min = 15.0
     
! xlv
      s1_bar    = s1_min + DTdz * z ! s1--> perturbation density [kg/m3]
      s1_bar_z  = dTdz            ! d/dz s1_bar  [kg/m4]
      s1_bar_zz = 0.0               ! dzz s1_bar  [kg/m5]
 

      return
      end subroutine eval_s1_bar
!=========================================================================


         subroutine eval_s2_bar(z,Lz,s2_bar,s2_bar_z,s2_bar_zz) 
!      input
!      z            vertical z spatial position  [m]
!      Lz           vertical size of computational domain [m]
!           
!     output 
!      s2_bar, s2_bar_z, s2_bar_zz
!      scalar concentration and z derivs   [dimensional units]

      implicit none
      real(kind=8) s2_bar,s2_bar_z,s2_bar_zz,z,Lz
!     LOCAL variables: you may need some intermediate variables,
!                      these need to be explicitly declared here
      real(kind=8)   :: s2_max=35.d0
 

          s2_bar = s2_max
          s2_bar_z = 0.0
          s2_bar_zz = 0.0       !     dzz s2_bar   [1/m2]

      return
     end subroutine eval_s2_bar



end module user_module