1 !=======================================================================
   2 !  global 1 deg model with 115 levels
   3 !======================================================================= 
   4 
   5 module config_module
   6  ! use this module only locally in this file
   7  implicit none
   8  real*8, allocatable :: t_star(:,:,:)
   9  real*8, allocatable :: s_star(:,:,:)
  10  real*8, allocatable :: qnec(:,:,:)
  11  real*8, allocatable :: qnet(:,:,:)
  12  real*8, allocatable :: qsol(:,:,:)
  13  real*8, allocatable :: divpen_shortwave(:)
  14  real*8, allocatable :: taux(:,:,:)
  15  real*8, allocatable :: tauy(:,:,:)
  16 end module config_module
  17 
  18 
  19 subroutine set_parameter
  20  ! ----------------------------------
  21  !       set here main parameter
  22  ! ----------------------------------
  23  use main_module   
  24  use config_module
  25  use eke_module   
  26  use tke_module   
  27  use idemix_module   
  28  use isoneutral_module   
  29  use diagnostics_module   
  30  implicit none
  31   nx   = 360
  32   ny   = 160
  33   nz   = 115
  34   dt_mom    = 3600.0!/2.0
  35   dt_tracer = 3600.0!/2.0
  36 
  37   coord_degree     = .true.
  38   enable_cyclic_x  = .true.
  39 
  40   runlen = 365.*86400*2
  41 
  42   enable_diag_ts_monitor = .true.; ts_monint = 86400.0
  43   enable_diag_snapshots  = .true.; snapint  = 365*86400.0!/12.
  44 
  45   enable_diag_overturning= .true.; overint = 365*86400; 
  46   overfreq = overint/24.
  47   enable_diag_energy     = .true.; energint = 365*86400; energfreq =overfreq
  48   enable_diag_averages   = .true.; aveint  = 365*86400; avefreq = overfreq
  49 
  50 
  51   congr_epsilon = 1e-6
  52   congr_max_iterations = 10000
  53   enable_streamfunction = .true.
  54   !enable_congrad_verbose = .true.
  55 
  56   enable_neutral_diffusion = .true.; 
  57       K_iso_0 = 1000.0
  58       K_iso_steep = 50.0
  59       iso_dslope=0.005
  60       iso_slopec=0.005
  61   enable_skew_diffusion = .true.
  62 
  63   enable_hor_friction = .true.; A_h = 5e4; 
  64   enable_hor_friction_cos_scaling = .true.; 
  65   hor_friction_cosPower=1
  66   enable_tempsalt_sources = .true.
  67 
  68   enable_implicit_vert_friction = .true.
  69   enable_tke = .true.
  70   c_k = 0.1
  71   c_eps = 0.7
  72   alpha_tke = 30.0
  73   mxl_min = 1d-8
  74   tke_mxl_choice = 2
  75   enable_tke_superbee_advection = .true.
  76 
  77   !K_gm_0 = 1000.0
  78   enable_eke = .true.
  79   eke_k_max  = 1e4
  80   eke_c_k    = 0.4
  81   eke_c_eps  = 0.5
  82   eke_cross  = 2.
  83   eke_crhin  = 1.0
  84   eke_lmin   = 100.0
  85   enable_eke_superbee_advection = .true.
  86   enable_eke_isopycnal_diffusion = .true.
  87 
  88   enable_idemix = .true.
  89   enable_eke_diss_surfbot = .true.
  90   eke_diss_surfbot_frac = 0.2 ! fraction which goes into bottom
  91   enable_idemix_superbee_advection = .true.
  92 
  93   enable_idemix_hor_diffusion = .true.; 
  94   !np=17+2
  95   !enable_idemix_M2 = .true.
  96   !enable_idemix_niw = .true.
  97   !omega_M2 =  2*pi/( 12*60*60 +  25.2 *60 )   ! M2 frequency in 1/s
  98 
  99   eq_of_state_type = 5
 100 end subroutine set_parameter
 101 
 102 
 103 
 104 subroutine set_grid
 105  use main_module   
 106  use config_module   
 107  implicit none
 108  real*4 :: dz4(nz)
 109  open(10,file='dz.bin',access='direct',recl=4*nz,status='old')
 110  read(10,rec=1) dz4
 111  close(10);
 112  dzt = dz4(nz:1:-1)
 113  dxt = 1.0
 114  dyt = 1.0
 115  y_Origin=-79.
 116  x_Origin=91.
 117 
 118 end subroutine set_grid
 119 
 120 
 121 subroutine set_coriolis
 122  use main_module   
 123  use config_module   
 124  implicit none
 125  integer :: j
 126  do j=js_pe-onx,je_pe+onx
 127    coriolis_t(:,j) = 2*omega*sin( yt(j)/180.*pi ) 
 128  enddo
 129 end subroutine set_coriolis
 130 
 131 
 132 subroutine set_initial_conditions
 133  use main_module   
 134  use config_module   
 135  use eke_module   
 136  use tke_module   
 137  use idemix_module   
 138  use isoneutral_module   
 139  implicit none
 140  integer :: i,j,k,kk,n
 141  real*4 :: dat4(nx,ny)
 142  real*8 :: pen(0:nz),swarg1,swarg2
 143  real*8 ::  rpart_shortwave  = 0.58
 144  real*8 ::  efold1_shortwave = 0.35
 145  real*8 ::  efold2_shortwave = 23.0
 146  include "netcdf.inc"
 147  integer :: iret, ncid,id
 148 
 149  allocate( t_star(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); t_star=0.0
 150  allocate( s_star(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); s_star=0.0
 151  allocate( qnec(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); qnec=0.0
 152  allocate( qnet(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); qnet=0.0
 153  allocate( qsol(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); qsol=0.0
 154  allocate( divpen_shortwave(nz) ); divpen_shortwave=0.0
 155  allocate( taux(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); taux=0.0
 156  allocate( tauy(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); tauy=0.0
 157 
 158  ! initital conditions
 159  open(10,file='lev_clim_temp.bin',access='direct',recl=4*nx*ny,status='old')
 160  open(20,file='lev_clim_salt.bin',access='direct',recl=4*nx*ny,status='old')
 161  kk=nz
 162  do k=1,nz
 163   read(10,rec=kk) dat4
 164   temp(is_pe:ie_pe,js_pe:je_pe,k,taum1) = dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,k)
 165   temp(is_pe:ie_pe,js_pe:je_pe,k,tau)   = dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,k)
 166   read(20,rec=kk) dat4
 167   salt(is_pe:ie_pe,js_pe:je_pe,k,tau)   = dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,k)
 168   salt(is_pe:ie_pe,js_pe:je_pe,k,taum1) = dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,k)
 169   kk=kk-1
 170  enddo
 171  close(10); close(20)
 172 
 173  ! wind stress on MIT grid
 174  open(10,file='ECMWFBB_taux.bin',access='direct',recl=4*nx*ny,status='old')
 175  open(20,file='ECMWFBB_tauy.bin',access='direct',recl=4*nx*ny,status='old')
 176  do n=1,12
 177   read(10,rec=n) dat4
 178   taux(is_pe:ie_pe,js_pe:je_pe,n)= dat4(is_pe:ie_pe,js_pe:je_pe)/rho_0
 179   read(20,rec=n) dat4
 180   tauy(is_pe:ie_pe,js_pe:je_pe,n)= dat4(is_pe:ie_pe,js_pe:je_pe)/rho_0
 181   call border_exchg_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,taux(:,:,n)) 
 182   call setcyclic_xy   (is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,taux(:,:,n))
 183   call border_exchg_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,tauy(:,:,n)) 
 184   call setcyclic_xy   (is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,tauy(:,:,n))
 185  enddo
 186  close(10); close(20)
 187  ! check for special values
 188  where( taux < -99.9) taux = 0
 189  where( tauy < -99.9) tauy = 0
 190 
 191 
 192  ! Qnet and dQ/dT and Qsol
 193  open(10,file='ECMWFBB_qnet.bin',access='direct',recl=4*nx*ny,status='old')
 194  open(20,file='ECMWFBB_dqdt.bin',access='direct',recl=4*nx*ny,status='old')
 195  open(30,file='ECMWFBB_swf.bin', access='direct',recl=4*nx*ny,status='old')
 196  do n=1,12
 197   read(10,rec=n) dat4
 198   qnet(is_pe:ie_pe,js_pe:je_pe,n) = -dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,nz)
 199   read(20,rec=n) dat4
 200   qnec(is_pe:ie_pe,js_pe:je_pe,n) = +dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,nz)
 201   read(30,rec=n) dat4
 202   qsol(is_pe:ie_pe,js_pe:je_pe,n) = -dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,nz)
 203  enddo
 204  close(10); close(20); close(30)
 205 
 206 
 207  ! SST and SSS
 208  open(10,file='ECMWFBB_target_sst.bin',access='direct',recl=4*nx*ny,status='old')
 209  open(20,file='lev_sss.bin',access='direct',recl=4*nx*ny,status='old')
 210  do n=1,12
 211   read(10,rec=n) dat4
 212   t_star(is_pe:ie_pe,js_pe:je_pe,n) = dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,nz)
 213   read(20,rec=n) dat4
 214   s_star(is_pe:ie_pe,js_pe:je_pe,n) = dat4(is_pe:ie_pe,js_pe:je_pe)*maskT(is_pe:ie_pe,js_pe:je_pe,nz)
 215  enddo
 216  close(10); close(20)
 217 
 218 
 219  if (enable_idemix ) then
 220 
 221 
 222   open(10,file='tidal_energy.bin',access='direct',recl=4*nx*ny,status='old')
 223   read(10,rec=1) dat4
 224   close(10)
 225   do j=js_pe,je_pe
 226    do i=is_pe,ie_pe
 227       k=max(1,kbot(i,j))
 228       dat4(i,j) = dat4(i,j)*maskW(i,j,k)/rho_0
 229    enddo
 230   enddo
 231 
 232   if (enable_idemix_M2) then
 233    do j=js_pe,je_pe
 234      do k=2,np-1
 235       forc_M2(is_pe:ie_pe,j,k) = 0.5*dat4(is_pe:ie_pe,j)/(2*pi)
 236      enddo
 237      forc_iw_bottom(is_pe:ie_pe,j) = 0.5*dat4(is_pe:ie_pe,j)
 238    enddo
 239   else
 240    forc_iw_bottom(is_pe:ie_pe,js_pe:je_pe) = dat4(is_pe:ie_pe,js_pe:je_pe)
 241   endif
 242 
 243   open(10,file='wind_energy.bin',access='direct',recl=4*nx*ny,status='old')
 244   read(10,rec=1) dat4
 245   close(10)
 246   do j=js_pe,je_pe
 247    do i=is_pe,ie_pe
 248       dat4(i,j) = dat4(i,j)*maskW(i,j,nz)/rho_0*0.2
 249    enddo
 250   enddo
 251 
 252   if (enable_idemix_niw) then
 253    do j=js_pe,je_pe
 254      do k=2,np-1
 255       forc_niw(is_pe:ie_pe,j,k) = 1.0*dat4(is_pe:ie_pe,j)/(2*pi)
 256      enddo
 257      forc_iw_surface(is_pe:ie_pe,j) = 0.0*dat4(is_pe:ie_pe,j)
 258    enddo
 259   else
 260    forc_iw_surface(is_pe:ie_pe,js_pe:je_pe) = dat4(is_pe:ie_pe,js_pe:je_pe)
 261   endif
 262 
 263   if (enable_idemix_niw) then
 264     do j=js_pe-onx,je_pe+onx
 265      omega_niw(j) = max(1d-8, abs( 1.05 * coriolis_t(j) ) )
 266     enddo
 267   endif
 268 
 269   if (enable_idemix_niw .or. enable_idemix_M2) then
 270    iret = nf_open('hrms_1deg.nc',NF_nowrite,ncid)
 271    iret = nf_inq_varid(ncid,'HRMS',id)
 272    iret = nf_get_vara_double(ncid,id ,(/is_pe,js_pe/),(/ie_pe-is_pe+1,je_pe-js_pe+1/),topo_hrms(is_pe:ie_pe,js_pe:je_pe))
 273    iret = nf_inq_varid(ncid,'LAM',id)
 274    iret = nf_get_vara_double(ncid,id ,(/is_pe,js_pe/),(/ie_pe-is_pe+1,je_pe-js_pe+1/),topo_lam(is_pe:ie_pe,js_pe:je_pe))
 275    call ncclos (ncid, iret)
 276 
 277    call border_exchg_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,topo_hrms) 
 278    call setcyclic_xy   (is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,topo_hrms)
 279    call border_exchg_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,topo_lam) 
 280    call setcyclic_xy   (is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,topo_lam)
 281   endif
 282 
 283  endif
 284 
 285 !     Initialize penetration profile for solar radiation
 286 !     and store divergence in divpen
 287 !     note that pen(nz) is set 0.0 instead of 1.0 to compensate for the
 288 !     shortwave part of the total surface flux 
 289   pen = 0.0
 290   do k=1,nz-1
 291      swarg1 = zw(k)/efold1_shortwave
 292      swarg2 = zw(k)/efold2_shortwave
 293      pen(k) = rpart_shortwave*exp(swarg1)+ (1.0-rpart_shortwave)*exp(swarg2)
 294   enddo
 295   do k=1,nz
 296      divpen_shortwave(k) = (pen(k) - pen(k-1))/dzt(k)
 297   enddo
 298 end subroutine set_initial_conditions
 299 
 300 
 301 
 302 subroutine get_periodic_interval(currentTime,cycleLength,recSpacing,nbrec,trec1,trec2,wght1,wght2)
 303  ! interpolation routine taken from mitgcm
 304  implicit none
 305  real*8, intent(in) :: currentTime,recSpacing,cycleLength
 306  integer, intent(in) :: nbrec
 307  real*8, intent(out) :: wght1,wght2
 308  integer, intent(out) :: tRec1,tRec2
 309  real*8 :: locTime,tmpTime
 310  locTime = currentTime - recSpacing*0.5 + cycleLength*( 2 - NINT(currentTime/cycleLength) )
 311  tmpTime = MOD( locTime, cycleLength )
 312  tRec1 = 1 + INT( tmpTime/recSpacing )
 313  tRec2 = 1 + MOD( tRec1, nbRec )
 314  wght2 = ( tmpTime - recSpacing*(tRec1 - 1) )/recSpacing
 315  wght1 = 1d0 - wght2
 316 end subroutine
 317 
 318 
 319 
 320 subroutine set_forcing
 321  use main_module   
 322  use config_module   
 323  use tke_module   
 324  implicit none
 325  integer :: i,j,k,n1,n2
 326  real*8 :: t_rest= 30*86400, cp_0 = 3991.86795711963  ! J/kg /K
 327  real*8 :: f1,f2,fxa,qqnet,qqnec,ice(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx)
 328 
 329  fxa = 365*86400d0
 330  call get_periodic_interval((itt-1)*dt_tracer,fxa,fxa/12.,12,n1,n2,f1,f2)
 331 
 332  ! linearly interpolate wind stress and shift from MITgcm U/V grid to this grid
 333  do j=js_pe-onx,je_pe+onx-1
 334   do i=is_pe-onx,ie_pe+onx-1
 335    surface_taux(i,j) = f1*taux(i+1,j,n1) + f2*taux(i+1,j,n2)
 336    surface_tauy(i,j) = f1*tauy(i,j+1,n1) + f2*tauy(i,j+1,n2)
 337   enddo
 338  enddo
 339 
 340  if (enable_tke ) then
 341   do j=js_pe-onx+1,je_pe+onx
 342    do i=is_pe-onx+1,ie_pe+onx
 343      forc_tke_surface(i,j) = sqrt( (0.5*(surface_taux(i,j)+surface_taux(i-1,j)))**2  &
 344                                   +(0.5*(surface_tauy(i,j)+surface_tauy(i,j-1)))**2 )**(3./2.) 
 345    enddo
 346   enddo
 347  endif
 348 
 349  !  W/m^2 K kg/J m^3/kg = K m/s 
 350  do j=js_pe,je_pe 
 351   do i=is_pe,ie_pe
 352    fxa   =  f1*t_star(i,j,n1)+f2*t_star(i,j,n2)
 353    qqnec =  f1*qnec(i,j,n1)+f2*qnec(i,j,n2)
 354    qqnet =  f1*qnet(i,j,n1)+f2*qnet(i,j,n2)
 355    forc_temp_surface(i,j)=(qqnet+qqnec*(fxa -temp(i,j,nz,tau)) )*maskT(i,j,nz)/cp_0/rho_0
 356 
 357    fxa  =  f1*s_star(i,j,n1)+f2*s_star(i,j,n2)
 358    forc_salt_surface(i,j)=dzt(nz)/t_rest*(fxa-salt(i,j,nz,tau))*maskT(i,j,nz)
 359   enddo
 360  enddo
 361 
 362  ! apply simple ice mask  
 363  ice = 1.0
 364  do j=js_pe,je_pe 
 365   do i=is_pe,ie_pe
 366    if (temp(i,j,nz,tau)*maskT(i,j,nz) <= -1.8 .and. forc_temp_surface(i,j) <=0 ) then
 367        forc_temp_surface(i,j) = 0.0
 368        forc_salt_surface(i,j) = 0.0
 369        ice(i,j)               = 0.0
 370    endif
 371   enddo
 372  enddo
 373 
 374  ! solar radiation
 375  do k=1,nz
 376   do j=js_pe,je_pe 
 377    do i=is_pe,ie_pe
 378     temp_source(i,j,k) =  (f1*qsol(i,j,n1)+f2*qsol(i,j,n2))*divpen_shortwave(k)*ice(i,j)*maskT(i,j,k)/cp_0/rho_0
 379    enddo
 380   enddo
 381  enddo
 382 end subroutine set_forcing
 383 
 384 
 385 
 386 subroutine set_topography
 387  use main_module   
 388  implicit none
 389  integer :: i,j,k,kk
 390  real*4 :: lev(nx,ny,nz),bathy(nx,ny)
 391  kbot=0
 392  open(10,file='bathymetry.bin',access='direct',recl=4*nx*ny,status='old')
 393  read(10,rec=1) bathy
 394  close(10)
 395 
 396  open(10,file='lev_clim_salt.bin',access='direct',recl=4*nx*ny,status='old')
 397  kk=nz
 398  do k=1,nz
 399   read(10,rec=kk) lev(:,:,k)
 400   kk=kk-1
 401  enddo
 402  close(10)
 403 
 404  do j=js_pe,je_pe
 405   do i=is_pe,ie_pe
 406    do k=nz,1,-1
 407     if (lev(i,j,k) .ne. 0.0 ) kbot(i,j)=k
 408    enddo
 409   enddo
 410   if (bathy(i,j) == 0.0) kbot(i,j) =0
 411  enddo
 412  where ( kbot == nz) kbot = 0
 413  !where (bathy == 0.0) kbot(1:nx,1:ny) =0
 414 
 415  do i=208,214
 416   do j=1,5
 417     if (i>=is_pe.and.i<=ie_pe.and.j>=js_pe.and.j<=je_pe) kbot(i,j)=0
 418   enddo
 419  enddo
 420  !kbot(208:214,1:5)=0
 421  i=105; j=135
 422  if (i>=is_pe.and.i<=ie_pe.and.j>=js_pe.and.j<=je_pe) kbot(i,j)=0
 423  !kbot(105,135)=0 ! Aleuten island
 424  do i=270,271
 425   j=131
 426     if (i>=is_pe.and.i<=ie_pe.and.j>=js_pe.and.j<=je_pe) kbot(i,j)=0
 427  enddo
 428  !kbot(270:271,131) = 0 ! Engl Kanal
 429 end subroutine set_topography
 430 
 431 
 432 
 433 
 434 
 435 
 436 
 437 subroutine set_diagnostics
 438  use main_module   
 439  use eke_module   
 440  use tke_module   
 441  use idemix_module   
 442  use isoneutral_module   
 443  implicit none
 444  call register_average('taux','Zonal wind stress','m^2/s','UT',surface_taux,0D0,.false.)
 445  call register_average('tauy','Meridional wind stress','m^2/s','TU',surface_tauy,0D0,.false.)
 446  call register_average('forc_temp_surface','Surface temperature flux','m K/s','TT',forc_temp_surface,0D0,.false.)
 447  call register_average('forc_salt_surface','Surface salinity flux','m g/s kg','TT',forc_salt_surface,0D0,.false.)
 448  if (enable_streamfunction) then
 449    call register_average('psi','Barotropic streamfunction','m^2/s','UU',psi(:,:,tau),0D0,.false.)
 450  else
 451    call register_average('psi','Surface pressure','m^2/s','TT',psi(:,:,tau),0D0,.false.)
 452  endif
 453  call register_average('temp','Temperature','deg C','TTT',0d0,temp(:,:,:,tau),.true.)
 454  call register_average('salt','Salinity','g/kg','TTT',0d0,salt(:,:,:,tau),.true.)
 455  call register_average('u','Zonal velocity','m/s','UTT',0d0,u(:,:,:,tau),.true.)
 456  call register_average('v','Meridional velocity','m/s','TUT',0d0,v(:,:,:,tau),.true.)
 457  call register_average('w','Vertical velocity','m/s','TTU',0d0,w(:,:,:,tau),.true.)
 458  call register_average('Nsqr','Square of stability frequency','1/s^2','TTU',0d0,Nsqr(:,:,:,tau),.true.)
 459  call register_average('Hd','Dynamic enthalpy','m^2/s^2','TTT',0d0,Hd(:,:,:,tau),.true.)
 460  call register_average('rho','Density','kg/m^3','TTT',0d0,rho(:,:,:,tau),.true.)
 461 
 462  call register_average('K_diss_v','Dissipation by vertical friction','m^2/s^3','TTU',0d0,K_diss_v,.true.)
 463  call register_average('K_diss_h','Dissipation by lateral friction','m^2/s^3','TTU',0d0,K_diss_h,.true.)
 464  call register_average('K_diss_bot','Dissipation by bottom friction','m^2/s^3','TTU',0d0,K_diss_bot,.true.)
 465  call register_average('P_diss_v','Dissipation by vertical mixing','m^2/s^3','TTU',0d0,P_diss_v,.true.)
 466  call register_average('P_diss_nonlin','Dissipation by nonlinear vert. mix.','m^2/s^3','TTU',0d0,P_diss_nonlin,.true.)
 467  call register_average('P_diss_iso','Dissipation by Redi mixing tensor','m^2/s^3','TTU',0d0,P_diss_iso,.true.)
 468 
 469  call register_average('kappaH','Vertical diffusivity','m^2/s','TTU',0d0,kappaH,.true.)
 470  if (enable_skew_diffusion)  then
 471    call register_average('B1_gm','Zonal component of GM streamfct.','m^2/s','TUT',0d0,B1_gm,.true.)
 472    call register_average('B2_gm','Meridional component of GM streamfct.','m^2/s','UTT',0d0,B2_gm,.true.)
 473  endif
 474  if (enable_TEM_friction)  then
 475    call register_average('kappa_gm','Vertical GM viscosity','m^2/s','TTU',0d0,kappa_gm,.true.)
 476    call register_average('K_diss_gm','Dissipation by GM friction','m^2/s^3','TTU',0d0,K_diss_gm,.true.)
 477  endif
 478 
 479  if (enable_tke)  then
 480    call register_average('TKE','Turbulent kinetic energy','m^2/s^2','TTU',0d0,tke(:,:,:,tau),.true.)
 481    call register_average('Prandtl','Prandtl number',' ','TTU',0d0,Prandtlnumber,.true.)
 482    call register_average('mxl','Mixing length',' ','TTU',0d0,mxl,.true.)
 483    call register_average('tke_diss','Dissipation of TKE','m^2/s^3','TTU',0d0,tke_diss,.true.)
 484    call register_average('forc_tke_surface','TKE surface forcing','m^3/s^2','TT',forc_tke_surface,0D0,.false.)
 485    call register_average('tke_surface_corr','TKE surface flux correction','m^3/s^2','TT',tke_surf_corr,0D0,.false.)
 486  endif
 487  if (enable_idemix)  then
 488    call register_average('E_iw','Internal wave energy','m^2/s^2','TTU',0d0,e_iw(:,:,:,tau),.true.)
 489    call register_average('forc_iw_surface','IW surface forcing','m^3/s^2','TT',forc_iw_surface,0D0,.false.)
 490    call register_average('forc_iw_bottom','IW bottom forcing','m^3/s^2','TT',forc_iw_bottom,0D0,.false.)
 491    call register_average('iw_diss','Dissipation of E_iw','m^2/s^3','TTU',0d0,iw_diss,.true.)
 492    call register_average('c0','Vertical IW group velocity','m/s','TTU',0d0,c0,.true.)
 493    call register_average('v0','Horizontal IW group velocity','m/s','TTU',0d0,v0,.true.)
 494  endif
 495  if (enable_eke)  then
 496   call register_average('EKE','Eddy energy','m^2/s^2','TTU',0d0,eke(:,:,:,tau),.true.)
 497   call register_average('K_gm','Lateral diffusivity','m^2/s','TTU',0d0,K_gm,.true.)
 498   call register_average('eke_diss','Eddy energy dissipation','m^2/s^3','TTU',0d0,eke_diss,.true.)
 499   call register_average('L_Rossby','Rossby radius','m','TT',L_rossby,0d0,.false.)
 500   call register_average('L_Rhines','Rhines scale','m','TTU',0d0,L_Rhines,.true.)
 501  endif
 502 
 503  if (enable_idemix_M2)  then
 504    call register_average('E_M2','M2 tidal energy','m^2/s^2','TT',E_M2_int,0d0,.false.)
 505    call register_average('cg_M2','M2 group velocity','m/s','TT',cg_M2,0d0,.false.)
 506    call register_average('tau_M2','Decay scale','1/s','TT',tau_M2,0d0,.false.)
 507    call register_average('alpha_M2_cont','Interaction coeff.','s/m^3','TT',alpha_M2_cont,0d0,.false.)
 508  endif
 509 
 510  if (enable_idemix_niw)  then
 511    call register_average('E_niw','NIW energy','m^2/s^2','TT',E_niw_int,0d0,.false.)
 512    call register_average('cg_niw','NIW group velocity','m/s','TT',cg_niw,0d0,.false.)
 513    call register_average('tau_niw','Decay scale','1/s','TT',tau_niw,0d0,.false.)
 514  endif
 515 end subroutine set_diagnostics