1 !=======================================================================
   2 ! Wide eddying channel with restoring zones at side walls
   3 ! Experiment CHANNEL in Eden (2010), Ocean modeling 32 (2010) 58-71
   4 !=======================================================================
   5 
   6 module config_module
   7  ! use this module only locally in this file
   8  implicit none
   9  real*8,parameter :: hRESOLVE = 0.5  ! 2 in original model 
  10  real*8,parameter :: vRESOLVE = 0.5  ! 2 in original model
  11  real*8,parameter :: M_0     = sqrt(1e-5*0.1/1024.*9.81)
  12  real*8,parameter :: N_0     = 0.004
  13  integer :: spg_width=int(3*hRESOLVE)
  14  real*8 :: T_rest=1./(5.*86400)
  15  !real*8 :: DX  = 32094.1729769
  16  !real*8,parameter :: DX  = 30e3
  17  !real*8,parameter :: Lx  = DX*128
  18  real*8,parameter :: Lx  = 3800e3
  19  real*8,parameter :: H   = 1800.0
  20 end module config_module
  21 
  22 
  23 
  24 
  25 
  26 subroutine set_parameter
  27  ! ----------------------------------
  28  !       set here main parameter
  29  ! ----------------------------------
  30  use main_module
  31  use config_module
  32  use diagnostics_module   
  33  use tke_module   
  34  use idemix_module   
  35  implicit none
  36 
  37  nx = int(128*hRESOLVE) ; nz = int(18*vRESOLVE); ny = int(128*hRESOLVE)
  38  dt_tracer = 1800.0/hRESOLVE
  39  dt_mom    = 1800.0/hRESOLVE
  40 
  41  coord_degree           = .false.
  42  enable_cyclic_x        = .true.
  43  enable_hydrostatic     = .true.
  44  eq_of_state_type       = 1
  45      
  46  congr_epsilon = 1e-12
  47  congr_max_iterations = 5000
  48  enable_streamfunction = .true.
  49      
  50  enable_ray_friction      = .true.
  51  r_ray = 0.5e-7 
  52      
  53  !enable_hor_friction = .false.
  54  !A_h = 100/HRESOLVE**2  
  55 
  56  enable_biharmonic_mixing    = .true.  ! was enabled in original model
  57  K_hbi  = 1e11/hRESOLVE**4
  58  !enable_superbee_advection = .true.
  59 
  60  enable_biharmonic_friction  = .true.  ! was enabled in original model
  61  A_hbi  = 1e11/hRESOLVE**4
  62      
  63  enable_tempsalt_sources =  .true.
  64 
  65  enable_conserve_energy = .false.
  66  kappah_0=1e-4/VRESOLVE**2
  67  kappam_0=1e-3/VRESOLVE**2
  68 
  69  !enable_implicit_vert_friction = .true.;
  70  !enable_tke = .true.
  71  c_k = 0.1
  72  c_eps = 0.7
  73  alpha_tke = 30.0
  74  mxl_min = 1d-8
  75  tke_mxl_choice = 2
  76  enable_tke_superbee_advection = .true.
  77 
  78  !enable_idemix = .true.
  79  enable_idemix_hor_diffusion = .true.;
  80  enable_idemix_superbee_advection = .true.
  81 
  82  runlen =  365*86400.*10
  83  enable_diag_ts_monitor = .true.; ts_monint = dt_tracer
  84  enable_diag_snapshots  = .true.; snapint  =  3*86400
  85  !enable_diag_averages   = .true.
  86  !aveint  = 365*86400./12!*100
  87  !avefreq = dt_tracer*10
  88  !enable_diag_energy     = .true.; energint =   86400*3.
  89  !energfreq = dt_tracer*10
  90 
  91   enable_diag_particles = .true.; particles_int = 3*86400.0
  92 
  93 end subroutine set_parameter
  94 
  95 
  96 subroutine set_grid
  97  use main_module   
  98  use config_module   
  99  implicit none
 100  !dxt = 1./3.*degtom*cos(30./180.*pi)/hRESOLVE ! original grid
 101  !dzt = 100.0 /vRESOLVE
 102  dxt = Lx/nx
 103  dyt = Lx/ny
 104  dzt = H/nz
 105 end subroutine set_grid
 106 
 107 
 108 subroutine set_coriolis
 109  use main_module   
 110  use config_module
 111  implicit none
 112  integer :: j
 113  real*8 :: phi0 , betaloc
 114  phi0 = 10.0 /180. *pi
 115  betaloc = 2*omega*cos(phi0)/radius
 116  do j=js_pe-onx,je_pe+onx
 117     coriolis_t(:,j) = 2*omega*sin(phi0) +betaloc*yt(j)
 118  enddo
 119 end subroutine set_coriolis
 120 
 121 
 122 subroutine set_initial_conditions
 123  ! ----------------------------------
 124  !      add here initial conditions
 125  ! ----------------------------------
 126  use main_module
 127  use config_module
 128  implicit none
 129  integer :: i,j,k
 130  real*8 :: B0,alpha,get_drhodT
 131  do k=1,nz
 132    do j=js_pe,je_pe
 133      do i=is_pe,ie_pe
 134        alpha = get_drhodT(salt(i,j,k,tau),temp(i,j,k,tau),zt(k) ) 
 135        B0=M_0**2*yt(j)+0.5e-3*sin(xt(i)/Lx*8.5*pi)!*exp(-(y-0.5)**2/0.5**2)
 136        temp(i,j,k,:)  = ( 32+rho_0/grav/alpha*(B0-N_0**2*zt(k))  )*maskT(i,j,k)
 137      enddo
 138    enddo
 139  enddo
 140 end subroutine set_initial_conditions
 141 
 142 
 143 
 144 subroutine set_forcing
 145  ! ----------------------------------
 146  !      add here restoring zones
 147  ! ----------------------------------
 148  use main_module
 149  use config_module
 150  implicit none
 151  integer :: i,j,k 
 152  real*8 :: B0, alpha, get_drhodT
 153  if (enable_tempsalt_sources) then
 154   do k=1,nz
 155    do j=2,spg_width+1
 156      if (j>=js_pe .and.  j<=je_pe) then
 157        do i=is_pe,ie_pe
 158          alpha = get_drhodT(salt(i,j,k,tau),temp(i,j,k,tau),zt(k) ) 
 159          B0=32+( yt(j)*M_0**2-N_0**2*zt(k) )*rho_0/grav/alpha
 160          temp_source(i,j,k)=maskT(i,j,k)*t_rest/(j-1.)*(B0-temp(i,j,k,tau))
 161        enddo
 162       endif
 163    enddo
 164    do j=ny-1,ny-spg_width,-1
 165      if (j>=js_pe .and.  j<=je_pe) then
 166        do i=is_pe,ie_pe
 167          alpha = get_drhodT(salt(i,j,k,tau),temp(i,j,k,tau),zt(k) ) 
 168          B0=32+( yt(j)*M_0**2-N_0**2*zt(k) )*rho_0/grav/alpha
 169          temp_source(i,j,k)=maskT(i,j,k)*t_rest/(-1.*(j-ny))*(B0-temp(i,j,k,tau))
 170        enddo
 171       endif
 172    enddo
 173   enddo
 174  endif
 175 end subroutine set_forcing
 176 
 177 
 178 subroutine set_topography
 179  use main_module   
 180  implicit none
 181  kbot=1
 182 end subroutine set_topography
 183 
 184 
 185 subroutine set_diagnostics
 186  use main_module   
 187  use eke_module   
 188  use tke_module   
 189  use idemix_module   
 190  use isoneutral_module   
 191  implicit none
 192  call register_average('taux','Zonal wind stress','m^2/s','UT',surface_taux,0D0,.false.)
 193  call register_average('tauy','Meridional wind stress','m^2/s','TU',surface_tauy,0D0,.false.)
 194  call register_average('forc_temp_surface','Surface temperature flux','m K/s','TT',forc_temp_surface,0D0,.false.)
 195  call register_average('forc_salt_surface','Surface salinity flux','m g/s kg','TT',forc_salt_surface,0D0,.false.)
 196  if (enable_streamfunction) then
 197    call register_average('psi','Barotropic streamfunction','m^2/s','UU',psi(:,:,tau),0D0,.false.)
 198  else
 199    call register_average('psi','Surface pressure','m^2/s','TT',psi(:,:,tau),0D0,.false.)
 200  endif
 201  call register_average('temp','Temperature','deg C','TTT',0d0,temp(:,:,:,tau),.true.)
 202  call register_average('salt','Salinity','g/kg','TTT',0d0,salt(:,:,:,tau),.true.)
 203  call register_average('u','Zonal velocity','m/s','UTT',0d0,u(:,:,:,tau),.true.)
 204  call register_average('v','Meridional velocity','m/s','TUT',0d0,v(:,:,:,tau),.true.)
 205  call register_average('w','Vertical velocity','m/s','TTU',0d0,w(:,:,:,tau),.true.)
 206  call register_average('Nsqr','Square of stability frequency','1/s^2','TTU',0d0,Nsqr(:,:,:,tau),.true.)
 207  call register_average('Hd','Dynamic enthalpy','m^2/s^2','TTT',0d0,Hd(:,:,:,tau),.true.)
 208 
 209  call register_average('K_diss_v','Dissipation by vertical friction','m^2/s^3','TTU',0d0,K_diss_v,.true.)
 210  call register_average('K_diss_h','Dissipation by lateral friction','m^2/s^3','TTU',0d0,K_diss_h,.true.)
 211  call register_average('K_diss_bot','Dissipation by bottom friction','m^2/s^3','TTU',0d0,K_diss_bot,.true.)
 212  call register_average('P_diss_v','Dissipation by vertical mixing','m^2/s^3','TTU',0d0,P_diss_v,.true.)
 213  call register_average('P_diss_nonlin','Dissipation by nonlinear vert. mix.','m^2/s^3','TTU',0d0,P_diss_nonlin,.true.)
 214  call register_average('P_diss_iso','Dissipation by Redi mixing tensor','m^2/s^3','TTU',0d0,P_diss_iso,.true.)
 215 
 216  call register_average('kappaH','Vertical diffusivity','m^2/s','TTU',0d0,kappaH,.true.)
 217  if (enable_skew_diffusion)  then
 218    call register_average('B1_gm','Zonal component of GM streamfct.','m^2/s','TUT',0d0,B1_gm,.true.)
 219    call register_average('B2_gm','Meridional component of GM streamfct.','m^2/s','UTT',0d0,B2_gm,.true.)
 220  else
 221    call register_average('kappa_gm','Vertical GM viscosity','m^2/s','TTU',0d0,kappa_gm,.true.)
 222    call register_average('K_diss_gm','Dissipation by GM friction','m^2/s^3','TTU',0d0,K_diss_gm,.true.)
 223  endif
 224 
 225  if (enable_tke)  then
 226    call register_average('TKE','Turbulent kinetic energy','m^2/s^2','TTU',0d0,tke(:,:,:,tau),.true.)
 227    call register_average('Prandtl','Prandtl number',' ','TTU',0d0,Prandtlnumber,.true.)
 228    call register_average('mxl','Mixing length',' ','TTU',0d0,mxl,.true.)
 229    call register_average('tke_diss','Dissipation of TKE','m^2/s^3','TTU',0d0,tke_diss,.true.)
 230    call register_average('forc_tke_surface','TKE surface forcing','m^3/s^2','TT',forc_tke_surface,0D0,.false.)
 231    call register_average('tke_surface_corr','TKE surface flux correction','m^3/s^2','TT',tke_surf_corr,0D0,.false.)
 232  endif
 233  if (enable_idemix)  then
 234    call register_average('E_iw','Internal wave energy','m^2/s^2','TTU',0d0,e_iw(:,:,:,tau),.true.)
 235    call register_average('forc_iw_surface','IW surface forcing','m^3/s^2','TT',forc_iw_surface,0D0,.false.)
 236    call register_average('forc_iw_bottom','IW bottom forcing','m^3/s^2','TT',forc_iw_bottom,0D0,.false.)
 237    call register_average('iw_diss','Dissipation of E_iw','m^2/s^3','TTU',0d0,iw_diss,.true.)
 238    call register_average('c0','Vertical IW group velocity','m/s','TTU',0d0,c0,.true.)
 239    call register_average('v0','Horizontal IW group velocity','m/s','TTU',0d0,v0,.true.)
 240  endif
 241  if (enable_eke)  then
 242   call register_average('EKE','Eddy energy','m^2/s^2','TTU',0d0,eke(:,:,:,tau),.true.)
 243   call register_average('K_gm','Lateral diffusivity','m^2/s','TTU',0d0,K_gm,.true.)
 244   call register_average('eke_diss','Eddy energy dissipation','m^2/s^3','TTU',0d0,eke_diss,.true.)
 245   call register_average('L_Rossby','Rossby radius','m','TT',L_rossby,0d0,.false.)
 246   call register_average('L_Rhines','Rhines scale','m','TTU',0d0,L_Rhines,.true.)
 247 
 248  endif
 249 
 250 end subroutine set_diagnostics
 251 
 252 
 253 
 254 
 255 subroutine set_particles
 256  use main_module   
 257  use particles_module   
 258  implicit none
 259  integer :: n
 260  real :: fxa
 261  real*8 :: xs,xe,zs,ze,ys,ye
 262 
 263  call allocate_particles(1000)
 264  xs=0;xe=nx*dxt(is_pe);
 265  ys=dyt(js_pe);ye=(ny-2)*dyt(js_pe);
 266  !zs=zt(1);ze=zt(nz)
 267  zs=-500;ze=-500
 268  do n=1,nptraj
 269     call random_number(fxa)
 270     pxyz(1,n) = xs+fxa*(xe-xs)
 271     call random_number(fxa)
 272     pxyz(2,n) = ys+fxa*(ye-ys)
 273     call random_number(fxa)
 274     pxyz(3,n) = zs+fxa*(ze-zs)
 275  enddo
 276 end subroutine set_particles