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