source: trunk/00FlowSolve_PL/PROJECTS/NEW_SRC/ibr.f90 @ 4

subroutine ibr_1
#define DEBUG_LEVEL 1
!     Routine to calculate the Inertia, Buoyancy, Rotation terms of 
!     xi momentum equation, ignoring the pressure and dissipation terms
!
!     inputs:
!            U,v,W contravariant velocities (x,z) and horizontal velocity (y)
!            GM_uu the grid metric term in front of the U^2 term
!            GM_uw the grid metric term in front of the U*W term
!            GM_ww the grid metric term in front of the W^2 term
!            pd     perturbation density
!            Rot  the inverse Rossby number U/fL
!            Ri   the Richardson number (N2*L^2)/U^2
!            nx,ny,locnz  size of fields
!            dxi,deta,dzeta  regular grid spacings on the computational mesh
!            work    work array of at least (6*nx*ny*locnz,3*nx*ny*locnz) 
!                    for evaluation of DNS terms, no evaluation of DNS terms
!            x_zeta  derivative of cartesion x wrt grid coordinate zeta
!            y_eta   derivative of cartesian y wrt coordinate eta
!            z_zeta  derivative of cartesion z wrt grid coordinate zeta
!            det_J determinant of transformation matrix
!            BC_flag - derivative conditions and end values
!
!     output:
!            rhs1 an (nx,ny,locnz) array such that dU/dt = rhs - pressure term
!
 
 use mpi_parameters
 use grid_info
 use dependent_variables
 use pde_parameters
 use rhs_variables
 use intermediate_variables
 use boundary_information
 use counters_flags_etc
 use dimensional_scales
 use user_parameters

 implicit none
 include '../input/problem_size.h'

 integer                                :: flag
 integer                                :: ii,jj,kk
 real                                   :: end_values(2),h
 real, dimension(:,:,:,:), allocatable  :: work
 allocate ( work(nx,ny,locnz,3) )
  work(:,:,:,:) = 0.

#if DEBUG_LEVEL >= 1
 if(myid==0)write(0,*) 'hello world from ibr_1, Rot=',Rot
#endif

!Calculate  dU/dxi, dU/dy dU/dzeta and store in work arrays 1,2,3
 flag=BC_flag(1,1)
 h=Grid(0)%dxi
!!! write(0,*) 'call compact_ddx, myid = ', myid
 call compact_ddx(U,work(1,1,1,1),flag,end_values,nx,ny,locnz)

 flag=BC_flag(1,2)
 h=Grid(0)%deta
!!! write(0,*) 'call compact_ddy, myid = ', myid
 call compact_ddy(U,work(1,1,1,2),flag,end_values,nx,ny,locnz)

 do j=1,ny
  work(:,j,:,2)=work(:,j,:,2)/Grid(0)%y_eta(j)
 enddo

 flag=BC_flag(1,3)  
 h=Grid(0)%dzeta
!!! write(0,*) 'call compact_ddz, myid = ', myid
 call compact_ddz_mpi(U,work(1,1,1,3),flag,end_values,nx,ny,locnz)

! Add up the advective terms
!!!  write(0,*) 'at ADD 1, myid = ', myid
!!!  do kk=1,locnz
!!!  do jj=1,ny
!!!  do ii=1,nx
!!!    work(ii,jj,kk,1) = 0.
!!!    work(ii,jj,kk,2) = 0.
!!!    work(ii,jj,kk,3) = 0.
!!!    U(ii,jj,kk) = 0.
!!!    V(ii,jj,kk) = 0.
!!!    W(ii,jj,kk) = 0.
!!!    if(work(ii,jj,kk,3) .lt. 0.) work(ii,jj,kk,3) = 0.
!!!  if(ii.ge.257) then
!!!  write(0,*) 'MYID',myid,'  II,JJ,KK,WORK = ',ii,jj,kk,work(ii,jj,kk,1),work(ii,jj,kk,2),work(ii,jj,kk,3)
!!!  endif
!!!    rhs1(ii,jj,kk,MM0) = -(work(ii,jj,kk,1) * U(ii,jj,kk) + &
!!!                         work(ii,jj,kk,2) * V(ii,jj,kk) +   &
!!!                         work(ii,jj,kk,3) * W(ii,jj,kk) )
!!!  enddo
!!!  enddo
!!!  enddo
  rhs1(:,:,:,MM0) = -( work(:,:,:,1)*U + work(:,:,:,2)*V + work(:,:,:,3)*W )

! Add in the terms due to grid transformation
!!!  write(0,*) 'at ADD 2, myid = ', myid
  do j=1,ny
   rhs1(1:nx,j,1:locnz,MM0) =                                    &
     rhs1(1:nx,j,1:locnz,MM0) -                                  &
   ( GM_uu(1:nx,1:locnz,1)*U(1:nx,j,1:locnz)*U(1:nx,j,1:locnz)   &
   + GM_uw(1:nx,1:locnz,1)*U(1:nx,j,1:locnz)*W(1:nx,j,1:locnz)   &
   + GM_ww(1:nx,1:locnz,1)*W(1:nx,j,1:locnz)*W(1:nx,j,1:locnz) )
  enddo

! Add on the buoyancy and rotation terms
!!!  write(0,*) 'at ADD 3, myid = ', myid
!!  do j=1,ny
!!   rhs1(1:nx,j,1:locnz,MM0) =                                     &
!!   rhs1(1:nx,j,1:locnz,MM0)                                       &
!!+ Grid(0)%x_zeta(1:nx,1:locnz)*Ri*pd(1:nx,j,1:locnz)/Grid(0)%det_J(1:nx,1:locnz)  &
!!+ Grid(0)%z_zeta(1:nx,1:locnz)*Rot*v(1:nx,j,1:locnz)/Grid(0)%det_J(1:nx,1:locnz)
!!  enddo
  do j=1,ny
   rhs1(1:nx,j,1:locnz,MM0) =  rhs1(1:nx,j,1:locnz,MM0)  &
       + Ri * pd(1:nx,j,1:locnz) * Grid(0)%x_zeta(1:nx,1:locnz)/Grid(0)%det_J(1:nx,1:locnz) &
       + (Rot +  beta_window(j) ) * v(1:nx,j,1:locnz) * Grid(0)%z_zeta(1:nx,1:locnz)/Grid(0)%det_J(1:nx,1:locnz)
  enddo

  print*,'beta: max=',maxval(Rot+beta_window(:))

 deallocate( work )    
end subroutine ibr_1

subroutine ibr_2
#define DEBUG_LEVEL 1
!     Routine to calculate the Inertia, Buoyancy, Rotation terms of 
!     eta momentum equation, ignoring the pressure and dissipation terms
!
!     inputs:
!            U,v,W contravariant velocities (x,z) and horizontal velocity (y)
!            GM_uu the grid metric term in front of the U^2 term
!            GM_uw the grid metric term in front of the U*W term
!            GM_ww the grid metric term in front of the W^2 term
!            Rot  the inverse Rossby number U/fL
!            Ri   the Richardson number (N2*L^2)/U^2
!            nx,ny,locnz  size of fields
!            dxi,deta,dzeta  regular grid spacings on the computational mesh
!            work    work array of at least (6*nx*ny*locnz,3*nx*ny*locnz) 
!                    for evaluation of DNS terms, no evaluation of DNS terms
!            x_zeta  derivative of cartesion x wrt grid coordinate zeta
!            y_eta   derivative of cartesian y wrt coordinate eta
!            z_zeta  derivative of cartesion z wrt grid coordinate zeta
!            det_J determinant of transformation matrix
!            BC_flag - derivative conditions and end values
!
!     output:
!            rhs2 an (nx,ny,locnz) array such that dv/dt = rhs - pressure term
!
 
 use mpi_parameters
 use grid_info
 use dependent_variables
 use intermediate_variables
 use pde_parameters
 use user_parameters
 use rhs_variables
 use boundary_information
 use counters_flags_etc
 use dimensional_scales

 implicit none
 include '../input/problem_size.h'

 integer                                :: flag
 real                                   :: end_values(2),h
 real, dimension(:,:,:,:), allocatable  :: work
 allocate ( work(nx,ny,locnz,3) )
 work(:,:,:,:) = 0.
 
#if DEBUG_LEVEL >= 1
 if(myid==0) write(0,*) 'hello world from ibr_2'
#endif

!Calculate  dv/dxi, dv/dy dv/dzeta and store in work arrays 1,2,3
 flag=BC_flag(1,1)
 h=Grid(0)%dxi
 call compact_ddx(v,work(1,1,1,1),flag,end_values,nx,ny,locnz)

 flag=BC_flag(1,2)
 h=Grid(0)%deta
 call compact_ddy(v,work(1,1,1,2),flag,end_values,nx,ny,locnz)
 do j=1,ny
   work(:,j,:,2)=work(:,j,:,2)/Grid(0)%y_eta(j)
 enddo

 flag=BC_flag(1,3)  
 h=Grid(0)%dzeta
 call compact_ddz_mpi(v,work(1,1,1,3),flag,end_values,nx,ny,locnz)

! Add up the advective terms
  rhs2(:,:,:,MM0) = -( work(:,:,:,1)*U + work(:,:,:,2)*V + work(:,:,:,3)*W )
        
! Add on the buoyancy and rotation terms
!!! xlv ajout effet beta
!!  do j=1,ny
!!   rhs2(1:nx,j,1:locnz,MM0) = rhs2(1:nx,j,1:locnz,MM0)         &
!!   -Rot*( Grid(0)%x_xi(1:nx,1:locnz)*U(1:nx,j,1:locnz) +       &
!!          Grid(0)%x_zeta(1:nx,1:locnz)*W(1:nx,j,1:locnz) )
!!  enddo
  do j=1,ny
   rhs2(1:nx,j,1:locnz,MM0) = rhs2(1:nx,j,1:locnz,MM0)         &
   -( Rot + beta_window(j) ) *( Grid(0)%x_xi(1:nx,1:locnz)*U(1:nx,j,1:locnz)  &
   + Grid(0)%x_zeta(1:nx,1:locnz)*W(1:nx,j,1:locnz) )
  enddo

     
  print*,'beta: max=',maxval(Rot+beta_window(:))


 deallocate( work )    
end subroutine ibr_2

subroutine ibr_3
#define DEBUG_LEVEL 1
!     Routine to calculate the Inertia, Buoyancy, Rotation terms of 
!     zeta momentum equation, ignoring the pressure and dissipation terms
!
!     inputs:
!            U,v,W contravariant velocities (x,z) and horizontal velocity (y)
!            GM_uu the grid metric term in front of the U^2 term
!            GM_uw the grid metric term in front of the U*W term
!            GM_ww the grid metric term in front of the W^2 term
!            pd     perturbation density
!            Rot  the inverse Rossby number U/fL
!            Ri   the Richardson number (N2*L^2)/U^2
!            nx,ny,locnz  size of fields
!            dxi,deta,dzeta  regular grid spacings on the computational mesh
!            work    work array of at least (6*nx*ny*locnz,3*nx*ny*locnz) 
!                    for evaluation of DNS terms, no evaluation of DNS terms
!            x_zeta  derivative of cartesion x wrt grid coordinate zeta
!            y_eta   derivative of cartesian y wrt coordinate eta
!            z_zeta  derivative of cartesion z wrt grid coordinate zeta
!            det_J determinant of transformation matrix
!            BC_flag - derivative conditions and end values
!
!     output:
!            rhs3 an (nx,ny,locnz) array such that dU/dt = rhs - pressure term
!
 
 use mpi_parameters
 use grid_info
 use dependent_variables
 use intermediate_variables
 use pde_parameters
 use rhs_variables
 use boundary_information
 use counters_flags_etc
 use user_parameters
 use dimensional_scales

 implicit none
 include '../input/problem_size.h'

 integer                                :: flag
 real                                   :: end_values(2),h
 real, dimension(:,:,:,:), allocatable  :: work
 real, dimension(:), allocatable        :: mean_value
 allocate ( work(nx,ny,locnz,3) )
 work(:,:,:,:) = 0.
 allocate( mean_value(locnz) )
 mean_value(:) = 0.

#if DEBUG_LEVEL >= 1
 if(myid==0) write(0,*) 'hello world from ibr_3'
#endif

!Calculate  dW/dxi, dW/dy dW/dzeta and store in work arrays 1,2,3
 flag=BC_flag(1,1)
 h=Grid(0)%dxi
 call compact_ddx(W,work(1,1,1,1),flag,end_values,nx,ny,locnz)

 flag=BC_flag(1,2)
 h=Grid(0)%deta
 call compact_ddy(W,work(1,1,1,2),flag,end_values,nx,ny,locnz)
 do j=1,ny
   work(:,j,:,2)=work(:,j,:,2)/Grid(0)%y_eta(j)
 enddo

 flag=BC_flag(1,3)  
 h=Grid(0)%dzeta
 call compact_ddz_mpi(W,work(1,1,1,3),flag,end_values,nx,ny,locnz)

! Add up the advective terms
  rhs3(:,:,:,MM0) = -( work(:,:,:,1)*U + work(:,:,:,2)*V + work(:,:,:,3)*W )
        
! Add in the terms due to grid transformation
  do j=1,ny
   rhs3(1:nx,j,1:locnz,MM0) =                                    &
     rhs3(1:nx,j,1:locnz,MM0) -                                  &
   ( GM_uu(1:nx,1:locnz,2)*U(1:nx,j,1:locnz)*U(1:nx,j,1:locnz)   &
   + GM_uw(1:nx,1:locnz,2)*U(1:nx,j,1:locnz)*W(1:nx,j,1:locnz)   &
   + GM_ww(1:nx,1:locnz,2)*W(1:nx,j,1:locnz)*W(1:nx,j,1:locnz) )
  enddo

! Compute the xy mean of the buoyancy term
  do k=1,locnz
   mean_value(k) = 0.0
   do j=1,ny
    do i=1,nx
     mean_value(k) = mean_value(k) + Grid(0)%x_xi(i,k)*Ri*pd(i,j,k)/Grid(0)%det_J(i,k)
    enddo
   enddo
   mean_value(k) = mean_value(k)/float(nx*ny)
  enddo

! Add on the buoyancy (- the mean component) and rotation terms
!!! xlv ajout effet beta
!!  do j=1,ny
!!   do k=1,locnz
!!    rhs3(1:nx,j,k,MM0) =             &
!!      rhs3(1:nx,j,k,MM0)             &
!!      - ( Grid(0)%x_xi(1:nx,k)*Ri*pd(1:nx,j,k)/Grid(0)%det_J(1:nx,k) - mean_value(k) )     &
!!      -   Grid(0)%z_xi(1:nx,k)*Rot*v(1:nx,j,k)/Grid(0)%det_J(1:nx,k) 
!!   enddo
!!  enddo
  do j=1,ny
   do k=1,locnz
    rhs3(1:nx,j,k,MM0) = rhs3(1:nx,j,k,MM0)  &
      - ( Ri * Grid(0)%x_xi(1:nx,k)/Grid(0)%det_J(1:nx,k) * pd(1:nx,j,k) - mean_value(k) ) &
      - ( Rot + beta_window(j) )* v(1:nx,j,k)* Grid(0)%z_xi(1:nx,k)/Grid(0)%det_J(1:nx,k)
   enddo
  enddo

  print*,'beta: max=',maxval(Rot+beta_window(:))

 deallocate( work,mean_value )    
end subroutine ibr_3

subroutine ibr_4
!     Routine to calculate the rhs of the perturbation scalar 1 equation
!     ignoring the diffusive term.
!
!     inputs:
!            U,v,W contravariant velocities (x,z) and horizontal velocity (y)
!            s1     perturbation scalar 1, s1_bar; the ambient profile
!            nx,ny,locnz  size of fields
!            dxi,deta,dzeta  regular grid spacings on the computational mesh
!            x_zeta  derivative of cartesion x wrt grid coordinate zeta
!            y_eta   derivative of cartesian y wrt coordinate eta
!            z_zeta  derivative of cartesion z wrt grid coordinate zeta
!            det_J determinant of transformation matrix
!            BC_flag - derivative conditions and end values
!
!     output:
!            rhs4 an (nx,ny,locnz) array such that ds1/dt = rhs - diffusive term
!
 
 use mpi_parameters
 use grid_info
 use dependent_variables
 use pde_parameters
 use rhs_variables
 use boundary_information
 use counters_flags_etc
 use user_parameters

 implicit none
 include '../input/problem_size.h'

 integer                                :: flag
 real                                   :: end_values(2),h
 real, dimension(:,:,:,:), allocatable  :: work
 allocate (work(nx,ny,locnz,4))
 work(:,:,:,:) = 0.
 
#if DEBUG_LEVEL >= 1
 if(myid==0) write(0,*) 'hello world from ibr_4'
#endif

!Store total scalar in work(:,:,:,4)
 do j=1,ny
  work(1:nx,j,1:locnz,4)=s1(1:nx,j,1:locnz)+s1_bar(1:nx,1:locnz)
 enddo

!Calculate  ds1/dxi, ds1/dy ds1/dzeta and store in work arrays 1,2,3
 flag=BC_flag(1,1)
 h=Grid(0)%dxi
 call compact_ddx(work(:,:,:,4),work(1,1,1,1),flag,end_values,nx,ny,locnz)

 flag=BC_flag(1,2)
 h=Grid(0)%deta
 call compact_ddy(work(:,:,:,4),work(1,1,1,2),flag,end_values,nx,ny,locnz)
 do j=1,ny
   work(:,j,:,2)=work(:,j,:,2)/Grid(0)%y_eta(j)
 enddo

 flag=BC_flag(1,3)  
 h=Grid(0)%dzeta
 call compact_ddz_mpi(work(:,:,:,4),work(1,1,1,3),flag,end_values,nx,ny,locnz)

! Add up the advective terms
  rhs4(:,:,:,MM0) = -( work(:,:,:,1)*U + work(:,:,:,2)*V + work(:,:,:,3)*W )
        
 deallocate( work )    
end subroutine ibr_4

subroutine ibr_5
!     Routine to calculate the rhs of the perturbation scalar 2 equation
!     ignoring the diffusive term.
!
!     inputs:
!            U,v,W contravariant velocities (x,z) and horizontal velocity (y)
!            s2     perturbation scalar 1, s2_bar; the ambient profile
!            nx,ny,locnz  size of fields
!            dxi,deta,dzeta  regular grid spacings on the computational mesh
!            x_zeta  derivative of cartesion x wrt grid coordinate zeta
!            y_eta   derivative of cartesian y wrt coordinate eta
!            z_zeta  derivative of cartesion z wrt grid coordinate zeta
!            det_J determinant of transformation matrix
!            BC_flag - derivative conditions and end values
!
!     output:
!            rhs5 an (nx,ny,locnz) array such that ds1/dt = rhs - diffusive term
!
 
 use mpi_parameters
 use grid_info
 use dependent_variables
 use pde_parameters
 use rhs_variables
 use boundary_information
 use counters_flags_etc
 use user_parameters

 implicit none
 include '../input/problem_size.h'

 integer                                :: flag
 real                                   :: end_values(2),h
 real, dimension(:,:,:,:), allocatable  :: work
 allocate (work(nx,ny,locnz,4))
 work(:,:,:,:) = 0.
 
#if DEBUG_LEVEL >= 1
 if(myid==0) write(0,*) 'hello world from ibr_5'
#endif

!Store total scalar in work(:,:,:,4)
 do j=1,ny
  work(1:nx,j,1:locnz,4)=s2(1:nx,j,1:locnz)+s2_bar(1:nx,1:locnz)
 enddo

!Calculate  ds2/dxi, ds2/dy ds2/dzeta and store in work arrays 1,2,3
 flag=BC_flag(1,1)
 h=Grid(0)%dxi
 call compact_ddx(work(:,:,:,4),work(1,1,1,1),flag,end_values,nx,ny,locnz)

 flag=BC_flag(1,2)
 h=Grid(0)%deta
 call compact_ddy(work(:,:,:,4),work(1,1,1,2),flag,end_values,nx,ny,locnz)
 do j=1,ny
   work(:,j,:,2)=work(:,j,:,2)/Grid(0)%y_eta(j)
 enddo

 flag=BC_flag(1,3)  
 h=Grid(0)%dzeta
 call compact_ddz_mpi(work(:,:,:,4),work(1,1,1,3),flag,end_values,nx,ny,locnz)

! Add up the advective terms
  rhs5(:,:,:,MM0) = -( work(:,:,:,1)*U + work(:,:,:,2)*V + work(:,:,:,3)*W )
        
 deallocate( work )    
end subroutine ibr_5