1 !=======================================================================
   2 !  idealised global ocean, 2x2 deg and 15 vertical levels
   3 !======================================================================= 
   4 
   5 module config_module
   6  implicit none
   7  real*8 :: yt_start = -39.0
   8  real*8 :: yt_end   = 43
   9  real*8 :: yu_start = -40.0
  10  real*8 :: yu_end   = 42
  11 end module config_module
  12 
  13 
  14 subroutine set_parameter
  15  ! ----------------------------------
  16  !       set here main parameter
  17  ! ----------------------------------
  18  use main_module   
  19  use eke_module   
  20  use tke_module   
  21  use idemix_module   
  22  use isoneutral_module   
  23  use diagnostics_module   
  24  implicit none
  25   nx   = 30; nz   = 15; ny  = 42
  26   dt_mom    = 4800  
  27   dt_tracer = 86400/2.0
  28 
  29   coord_degree     = .true.
  30   enable_cyclic_x  = .true.
  31 
  32   runlen =  365*86400.*100
  33 
  34   enable_diag_ts_monitor = .true.; ts_monint = 365*86400./12.
  35   enable_diag_snapshots  = .true.; snapint  =  365*86400./12.0
  36   !enable_diag_overturning= .true.; overint  =  365*86400./48.0; overfreq = dt_tracer
  37   enable_diag_energy     = .true.; energint =   365*86400./48; energfreq = dt_tracer*10
  38   !enable_diag_averages   = .true.; aveint  = 365*86400.; avefreq = dt_tracer*10
  39   !enable_diag_particles = .true.; particles_int = snapint
  40 
  41   congr_epsilon = 1e-12
  42   congr_max_iterations = 5000
  43   enable_streamfunction = .true.
  44 
  45   !enable_hor_diffusion = .true.;  K_h=2000
  46   enable_neutral_diffusion = .true.; 
  47       K_iso_0 = 1000.0
  48       K_iso_steep = 500.0
  49       iso_dslope=0.005
  50       iso_slopec=0.01
  51   !enable_skew_diffusion = .true.
  52   enable_TEM_friction = .true.
  53 
  54   enable_hor_friction = .true.; A_h = (2*degtom)**3*2e-11    
  55   enable_hor_friction_cos_scaling = .true.; hor_friction_cosPower=1
  56 
  57   !enable_bottom_friction = .true.; r_bot = 1e-5
  58   !enable_bottom_friction_var = .true.; 
  59   enable_quadratic_bottom_friction = .true.; r_quad_bot = 2e-3
  60   enable_store_bottom_friction_tke = .true.; 
  61 
  62   enable_implicit_vert_friction = .true.; 
  63   enable_tke = .true.
  64   c_k = 0.1
  65   c_eps = 0.7
  66   alpha_tke = 30.0
  67   mxl_min = 1d-8
  68   tke_mxl_choice = 2
  69   enable_tke_superbee_advection = .false.
  70 
  71   K_gm_0 = 1000
  72   enable_eke = .true.
  73   eke_k_max  = 1e4
  74   eke_c_k    = 0.4
  75   eke_c_eps  = 0.5
  76   eke_cross  = 2.
  77   eke_crhin  = 1.0
  78   eke_lmin   = 100.0
  79   eke_r_bot = 2e-3
  80   enable_eke_superbee_advection = .true.
  81   enable_eke_isopycnal_diffusion = .true.
  82 
  83   enable_eke_leewave_dissipation = .true.
  84   eke_int_diss0 = 1./(10*86400.)
  85   c_lee0 = 5.
  86   eke_Ri0 = 300.
  87   eke_Ri1 = 50.
  88   alpha_eke = 20.0
  89 
  90 
  91 
  92   enable_idemix = .true.
  93   enable_idemix_hor_diffusion = .true.; 
  94   !enable_eke_diss_surfbot = .true.
  95   !eke_diss_surfbot_frac = 0.2
  96   enable_idemix_superbee_advection = .true.
  97   
  98   eq_of_state_type = 3 
  99 end subroutine set_parameter
 100 
 101 
 102 subroutine set_grid
 103   use main_module   
 104   implicit none
 105   real*8 :: ddz(15)  = (/50.,70.,100.,140.,190.,240.,290.,340.,390.,440.,490.,540.,590.,640.,690./)
 106   if (nz == 15) then
 107     dzt = ddz(15:1:-1)/2.5
 108   else
 109     call halt_stop('in set grid')
 110   endif
 111   dxt = 2.0
 112   dyt = 2.0
 113   x_origin=  0.0
 114   y_origin= -40.0
 115 end subroutine set_grid
 116 
 117 
 118 subroutine set_coriolis
 119  use main_module   
 120  implicit none
 121  integer :: j
 122  do j=js_pe-onx,je_pe+onx
 123    coriolis_t(:,j) = 2*omega*sin( yt(j)/180.*pi ) 
 124  enddo
 125 end subroutine set_coriolis
 126 
 127 
 128 subroutine set_initial_conditions
 129  use main_module   
 130  use eke_module   
 131  use tke_module   
 132  use idemix_module   
 133  use config_module   
 134  implicit none
 135  integer :: i,j,k
 136 
 137  do k=1,nz
 138   !salt(:,:,k,:) = 34+(1-zt(k)/zw(1))*2
 139   salt(:,:,k,:) = 35.0
 140   temp(:,:,k,:) =  (1-zt(k)/zw(1))*15
 141  enddo
 142 
 143  do j=js_pe-onx,je_pe+onx
 144   if ( yt(j)< -20.0) surface_taux(:,j) =.1e-3*sin(pi*(yu(j)-yu_start)/(-20.0-yt_start))*maskU(:,j,nz)
 145   if ( yt(j)> 10.0)  surface_taux(:,j) =.1e-3*(1-cos(2*pi*(yu(j)-10.0)/(yu_end-10.0)))*maskU(:,j,nz)
 146  enddo
 147 
 148  if (enable_tke ) then
 149   do j=js_pe-onx+1,je_pe+onx
 150    do i=is_pe-onx+1,ie_pe+onx
 151      forc_tke_surface(i,j) = sqrt( (0.5*(surface_taux(i,j)+surface_taux(i-1,j)))**2  &
 152                                   +(0.5*(surface_tauy(i,j)+surface_tauy(i,j-1)))**2 )**(3./2.) 
 153    enddo
 154   enddo
 155  endif
 156 
 157  if (enable_idemix ) then
 158    forc_iw_bottom =  1e-6
 159    forc_iw_surface = 0.1e-6
 160  endif
 161 
 162  if (enable_bottom_friction_var) then
 163   do j=js_pe-onx,je_pe+onx
 164      if ( yt(j)< -20.0) r_bot_var_u(:,j) = r_bot
 165      if ( yt(j)< -20.0) r_bot_var_v(:,j) = r_bot
 166   enddo
 167  endif
 168 
 169  if (enable_eke .and. enable_eke_leewave_dissipation ) then
 170   eke_topo_hrms = 100.0
 171   eke_topo_lam =1e3/0.2 !1e3/lam = 0.2 -> lam = 1e3/0.2
 172  endif
 173 
 174 end subroutine set_initial_conditions
 175 
 176 
 177 
 178 function t_star(j)
 179  use main_module   
 180  use config_module   
 181  implicit none
 182  integer :: j
 183  real*8 :: t_star
 184  t_star=15
 185  if (yt(j)<-20.0) t_star=15*(yt(j)-yt_start)/(-20.0-yt_start)
 186  if (yt(j)> 20.0) t_star=15*(1-(yt(j)-20)/(yt_end -20) )
 187 end function t_star
 188 
 189 
 190 function s_star(j)
 191  use main_module   
 192  use config_module   
 193  implicit none
 194  integer :: j
 195  real*8 :: s_star
 196  s_star=35
 197  if (yt(j)<-30.0) s_star=5*(yt(j)-yt_start)/(-30.0-yt_start)+30
 198  if (yt(j)> 30.0) s_star=5*(1-(yt(j)-30)/(yt_end-30) )+30
 199 end function s_star
 200 
 201 
 202 
 203 
 204 subroutine set_forcing
 205  use main_module   
 206  implicit none
 207  integer :: i,j
 208  real*8 :: t_rest,t_star,fxa,fxb!,s_star
 209  ! P-E = 50 mg/m^2/s  = 50 10^-6 kg/m^2/s , salt flux = - S (P-E)   kg/m^2 /s
 210  ! dS/dt (1/s) =  d/dz  F  , F (m/s ) = S (E-P) /rho ( m/s)
 211  real*8, parameter :: s_flux0 = -50e-6 *0.035
 212 
 213  t_rest=30*86400
 214  do j=js_pe-onx,je_pe+onx
 215     forc_temp_surface(:,j)=dzt(nz)/t_rest*(t_star(j)-temp(:,j,nz,tau)) 
 216     !forc_salt_surface(:,j)=dzt(nz)/t_rest*(s_star(j)-salt(:,j,nz,tau)) 
 217  enddo
 218 
 219 
 220  !forc_salt_surface = 0.0
 221  !do j=js_pe,je_pe
 222  ! if (yt(j)>30.0 .or. yt(j)< -30.0) forc_salt_surface(:,j)= s_flux0*maskT(:,j,nz)
 223  !enddo
 224  !fxa = 0.0; fxb = 0.0
 225  !do j=js_pe,je_pe
 226  ! do i=is_pe,ie_pe
 227  !   fxa = fxa + forc_salt_surface(i,j)*area_t(i,j)*maskT(i,j,nz)
 228  !   fxb = fxb + area_t(i,j)*maskT(i,j,nz)
 229  ! enddo
 230  !enddo
 231  !call global_sum(fxa); call global_sum(fxb)
 232  !fxa = fxa/fxb
 233  !forc_salt_surface = forc_salt_surface  -fxa
 234 
 235 end subroutine set_forcing
 236 
 237 
 238 
 239 subroutine set_topography
 240  use main_module   
 241  implicit none
 242  integer :: i,j
 243  
 244  kbot=0
 245  do i=is_pe,ie_pe
 246    do j=js_pe,je_pe
 247      if ( (yt(j)<-20.0).or.(xt(i)<59.0 .and. xt(i)>1.0 ))  kbot(i,j)=1
 248    enddo
 249  enddo
 250 end subroutine set_topography
 251 
 252 
 253 
 254 
 255 subroutine set_diagnostics
 256  use main_module   
 257  use eke_module   
 258  use tke_module   
 259  use idemix_module   
 260  use isoneutral_module   
 261  implicit none
 262  call register_average('taux','Zonal wind stress','m^2/s','UT',surface_taux,0D0,.false.)
 263  call register_average('tauy','Meridional wind stress','m^2/s','TU',surface_tauy,0D0,.false.)
 264  call register_average('forc_temp_surface','Surface temperature flux','m K/s','TT',forc_temp_surface,0D0,.false.)
 265  call register_average('forc_salt_surface','Surface salinity flux','m g/s kg','TT',forc_salt_surface,0D0,.false.)
 266  if (enable_streamfunction) then
 267    call register_average('psi','Barotropic streamfunction','m^2/s','UU',psi(:,:,tau),0D0,.false.)
 268  else
 269    call register_average('psi','Surface pressure','m^2/s','TT',psi(:,:,tau),0D0,.false.)
 270  endif
 271  call register_average('temp','Temperature','deg C','TTT',0d0,temp(:,:,:,tau),.true.)
 272  call register_average('salt','Salinity','g/kg','TTT',0d0,salt(:,:,:,tau),.true.)
 273  call register_average('u','Zonal velocity','m/s','UTT',0d0,u(:,:,:,tau),.true.)
 274  call register_average('v','Meridional velocity','m/s','TUT',0d0,v(:,:,:,tau),.true.)
 275  call register_average('w','Vertical velocity','m/s','TTU',0d0,w(:,:,:,tau),.true.)
 276  call register_average('Nsqr','Square of stability frequency','1/s^2','TTU',0d0,Nsqr(:,:,:,tau),.true.)
 277  call register_average('Hd','Dynamic enthalpy','m^2/s^2','TTT',0d0,Hd(:,:,:,tau),.true.)
 278 
 279  call register_average('K_diss_v','Dissipation by vertical friction','m^2/s^3','TTU',0d0,K_diss_v,.true.)
 280  call register_average('K_diss_h','Dissipation by lateral friction','m^2/s^3','TTU',0d0,K_diss_h,.true.)
 281  call register_average('K_diss_bot','Dissipation by bottom friction','m^2/s^3','TTU',0d0,K_diss_bot,.true.)
 282  call register_average('P_diss_v','Dissipation by vertical mixing','m^2/s^3','TTU',0d0,P_diss_v,.true.)
 283  call register_average('P_diss_nonlin','Dissipation by nonlinear vert. mix.','m^2/s^3','TTU',0d0,P_diss_nonlin,.true.)
 284  call register_average('P_diss_iso','Dissipation by Redi mixing tensor','m^2/s^3','TTU',0d0,P_diss_iso,.true.)
 285 
 286  call register_average('kappaH','Vertical diffusivity','m^2/s','TTU',0d0,kappaH,.true.)
 287  if (enable_skew_diffusion)  then
 288    call register_average('B1_gm','Zonal component of GM streamfct.','m^2/s','TUT',0d0,B1_gm,.true.)
 289    call register_average('B2_gm','Meridional component of GM streamfct.','m^2/s','UTT',0d0,B2_gm,.true.)
 290  else
 291    call register_average('kappa_gm','Vertical GM viscosity','m^2/s','TTU',0d0,kappa_gm,.true.)
 292    call register_average('K_diss_gm','Dissipation by GM friction','m^2/s^3','TTU',0d0,K_diss_gm,.true.)
 293  endif
 294 
 295  if (enable_tke)  then
 296    call register_average('TKE','Turbulent kinetic energy','m^2/s^2','TTU',0d0,tke(:,:,:,tau),.true.)
 297    call register_average('Prandtl','Prandtl number',' ','TTU',0d0,Prandtlnumber,.true.)
 298    call register_average('mxl','Mixing length',' ','TTU',0d0,mxl,.true.)
 299    call register_average('tke_diss','Dissipation of TKE','m^2/s^3','TTU',0d0,tke_diss,.true.)
 300    call register_average('forc_tke_surface','TKE surface forcing','m^3/s^2','TT',forc_tke_surface,0D0,.false.)
 301    call register_average('tke_surface_corr','TKE surface flux correction','m^3/s^2','TT',tke_surf_corr,0D0,.false.)
 302  endif
 303  if (enable_idemix)  then
 304    call register_average('E_iw','Internal wave energy','m^2/s^2','TTU',0d0,e_iw(:,:,:,tau),.true.)
 305    call register_average('forc_iw_surface','IW surface forcing','m^3/s^2','TT',forc_iw_surface,0D0,.false.)
 306    call register_average('forc_iw_bottom','IW bottom forcing','m^3/s^2','TT',forc_iw_bottom,0D0,.false.)
 307    call register_average('iw_diss','Dissipation of E_iw','m^2/s^3','TTU',0d0,iw_diss,.true.)
 308    call register_average('c0','Vertical IW group velocity','m/s','TTU',0d0,c0,.true.)
 309    call register_average('v0','Horizontal IW group velocity','m/s','TTU',0d0,v0,.true.)
 310  endif
 311  if (enable_eke)  then
 312   call register_average('EKE','Eddy energy','m^2/s^2','TTU',0d0,eke(:,:,:,tau),.true.)
 313   call register_average('K_gm','Lateral diffusivity','m^2/s','TTU',0d0,K_gm,.true.)
 314   !call register_average('eke_diss','Eddy energy dissipation','m^2/s^3','TTU',0d0,eke_diss,.true.)
 315   call register_average('L_Rossby','Rossby radius','m','TT',L_rossby,0d0,.false.)
 316   call register_average('L_Rhines','Rhines scale','m','TTU',0d0,L_Rhines,.true.)
 317 
 318  endif
 319 
 320 end subroutine set_diagnostics
 321 
 322 
 323 
 324 
 325 
 326 
 327 
 328 subroutine set_particles
 329  use main_module   
 330  use particles_module   
 331  implicit none
 332  integer :: n
 333  real :: fxa
 334  real*8 :: xs,xe,zs,ze,ys,ye
 335 
 336  call allocate_particles(1000)
 337  xs=20;xe=20;
 338  ys=-40;ye=40;
 339  zs=-200;ze=-200
 340  do n=1,nptraj
 341     call random_number(fxa)
 342     pxyz(1,n) = xs+fxa*(xe-xs)
 343     call random_number(fxa)
 344     pxyz(2,n) = ys+fxa*(ye-ys)
 345     call random_number(fxa)
 346     pxyz(3,n) = zs+fxa*(ze-zs)
 347  enddo
 348 end subroutine set_particles