1 !=======================================================================
   2 !  FLAME 4/3 deg model
   3 !======================================================================= 
   4 
   5 
   6 
   7 module config_module
   8  ! use this module only locally in this file
   9  implicit none
  10  real*8, allocatable :: t_clim(:,:,:)
  11  real*8, allocatable :: s_clim(:,:,:)
  12  real*8, allocatable :: t_rest(:,:,:)
  13  real*8, allocatable :: s_rest(:,:,:)
  14  real*8, allocatable :: taux(:,:,:),tauy(:,:,:)
  15  real*8, allocatable :: t_star(:,:,:,:),s_star(:,:,:,:),tscl(:,:,:)
  16 end module config_module
  17 
  18  
  19 
  20 subroutine set_parameter
  21  ! ----------------------------------
  22  !       set here main parameter
  23  ! ----------------------------------
  24  use main_module   
  25  use config_module
  26  use eke_module   
  27  use tke_module   
  28  use idemix_module   
  29  use isoneutral_module   
  30  use diagnostics_module   
  31  implicit none
  32   nx   = 87
  33   ny   = 89
  34   nz   = 45
  35   dt_mom    = 3600.0
  36   dt_tracer = 3600.0
  37 
  38   coord_degree     = .true.
  39   runlen = 365.*86400 *5
  40 
  41   enable_diag_ts_monitor = .true.; ts_monint = 86400.0
  42   enable_diag_snapshots  = .true.; snapint  =  365*86400. /20.0
  43   enable_diag_energy     = .true.; energint = dt_tracer !365*86400./12.0
  44   energfreq = dt_tracer
  45   !enable_diag_averages   = .true.
  46   !aveint  = 365*86400. /12.
  47   !avefreq = 86400.0
  48 
  49   congr_epsilon = 1e-6
  50   congr_max_iterations = 20000
  51   enable_streamfunction = .true.
  52   !enable_congrad_verbose = .true.
  53 
  54 !  enable_hor_diffusion = .true.; K_h = 2000.0
  55   enable_neutral_diffusion = .true.; 
  56       K_iso_0 = 1000.0
  57       K_iso_steep = 200.0
  58       iso_dslope=1./1000.0
  59       iso_slopec=4./1000.0
  60   enable_skew_diffusion = .true.
  61 
  62   enable_hor_friction = .true.; A_h = 5e4; 
  63   enable_hor_friction_cos_scaling = .true.; hor_friction_cosPower = 3
  64   enable_tempsalt_sources = .true.
  65 
  66   !kappaH_0 = 1e-4; kappaM_0 = 10e-4
  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 
  75   K_gm_0 = 1000.0
  76   !enable_eke = .true.
  77 
  78   enable_idemix = .true.
  79   enable_idemix_hor_diffusion = .true.; 
  80 
  81   eq_of_state_type = 5
  82 end subroutine set_parameter
  83 
  84 
  85 
  86 subroutine set_grid
  87  use main_module   
  88  use config_module   
  89  implicit none
  90  include "netcdf.inc"
  91  integer :: iret,id,ncid
  92 
  93  iret=nf_open('forcing.cdf',NF_NOWRITE,ncid)
  94  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
  95 
  96  !if (n_pes>1) then
  97  !  call halt_stop(' cannot read grid')
  98  !endif
  99  
 100  iret=nf_inq_varid(ncid,'dxtdeg',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 101  !iret= nf_get_vara_double(ncid,id,1, nx, dxt(1:nx) )
 102  !dxt(1-onx:0)=dxt(1); dxt(nx+1:nx+onx)=dxt(nx)
 103  iret= nf_get_vara_double(ncid,id,is_pe, ie_pe-is_pe+1, dxt(is_pe:ie_pe) )
 104 
 105  iret=nf_inq_varid(ncid,'dytdeg',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 106  !iret= nf_get_vara_double(ncid,id,1, ny, dyt(1:ny) )
 107  !dyt(1-onx:0)=dyt(1); dyt(ny+1:ny+onx)=dyt(ny)
 108  iret= nf_get_vara_double(ncid,id,js_pe, je_pe-js_pe+1, dyt(js_pe:je_pe) )
 109 
 110  iret=nf_inq_varid(ncid,'dzt',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 111  iret= nf_get_vara_double(ncid,id,1,nz, dzt )
 112  dzt(1:nz) = dzt(nz:1:-1)/100.0
 113 
 114  iret=nf_inq_varid(ncid,'xu',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 115  iret= nf_get_vara_double(ncid,id,1,1, x_origin )
 116 
 117  iret=nf_inq_varid(ncid,'yu',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 118  iret= nf_get_vara_double(ncid,id,1,1, y_origin )
 119 
 120  iret = nf_close (ncid)
 121 
 122 end subroutine set_grid
 123 
 124 
 125 subroutine set_coriolis
 126  use main_module   
 127  use config_module   
 128  implicit none
 129  integer :: j
 130  do j=js_pe-onx,je_pe+onx
 131    coriolis_t(:,j) = 2*omega*sin( yt(j)/180.*pi ) 
 132  enddo
 133 end subroutine set_coriolis
 134 
 135 
 136 
 137 
 138 
 139 subroutine set_topography
 140  use main_module   
 141  implicit none
 142  include "netcdf.inc"
 143  integer :: iret,id,ncid
 144  integer :: kmt(nx,ny),i,j
 145 
 146  iret=nf_open('forcing.cdf',NF_NOWRITE,ncid)
 147  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 148  iret=nf_inq_varid(ncid,'kmt',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 149  iret= nf_get_vara_int(ncid,id,(/1,1/),(/nx,ny/), kmt )
 150  iret = nf_close (ncid)
 151 
 152  kbot=1
 153  do j=js_pe,je_pe
 154    do i=is_pe,ie_pe
 155     kbot(i,j)=min(nz,nz-kmt(i,j)+1)
 156     if ( kmt(i,j) == 0) kbot(i,j)=0
 157    enddo
 158  enddo
 159 
 160 end subroutine set_topography
 161 
 162 
 163 
 164 
 165 
 166 subroutine set_initial_conditions
 167  use main_module   
 168  use config_module   
 169  use eke_module   
 170  use tke_module   
 171  use idemix_module   
 172  use isoneutral_module   
 173  implicit none
 174  include "netcdf.inc"
 175  integer :: iret,id,ncid
 176  integer :: i,j,k,n
 177  real*8 :: dat8(nx,ny,12)
 178  real*4 :: dat4(nx,ny)
 179 
 180  allocate( t_clim(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); t_clim=0.0
 181  allocate( s_clim(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); s_clim=0.0
 182  allocate( t_rest(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); t_rest=0.0
 183  allocate( s_rest(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); s_rest=0.0
 184  allocate( taux(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); taux=0.0
 185  allocate( tauy(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); tauy=0.0
 186  allocate( t_star(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,nz,12) ); t_star=0.0
 187  allocate( s_star(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,nz,12) ); s_star=0.0
 188  allocate( tscl(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,nz) ); tscl=0.0
 189 
 190  ! initial conditions
 191  iret=nf_open('forcing.cdf',NF_NOWRITE,ncid)
 192  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 193 
 194  iret=nf_inq_varid(ncid,'temp_ic',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 195  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/),(/ie_pe-is_pe+1,je_pe-js_pe+1,nz/), temp(is_pe:ie_pe,js_pe:je_pe,1:nz,1) )
 196  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 197 
 198  iret=nf_inq_varid(ncid,'salt_ic',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 199  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/),(/ie_pe-is_pe+1,je_pe-js_pe+1,nz/), salt(is_pe:ie_pe,js_pe:je_pe,1:nz,1) )
 200  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 201 
 202  temp(:,:,1:nz,1) = temp(:,:,nz:1:-1,1)
 203  salt(:,:,1:nz,1) = salt(:,:,nz:1:-1,1)
 204 
 205  do k=1,nz
 206    do j=js_pe,je_pe
 207     do i=is_pe,ie_pe
 208       temp(i,j,k,1:3) = temp(i,j,k,1)*maskT(i,j,k)
 209       salt(i,j,k,1:3) = (salt(i,j,k,1)*1000+35)*maskT(i,j,k)
 210     enddo
 211    enddo
 212  enddo
 213 
 214  ! wind stress
 215  iret=nf_inq_varid(ncid,'taux',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 216  iret=nf_get_vara_double(ncid,id,(/1,1,1/),(/nx,ny,12/), dat8 )
 217  where (dat8 <= -1e20 ) dat8 = 0
 218  do n=1,12
 219   taux(is_pe:ie_pe,js_pe:je_pe,n) = dat8(is_pe:ie_pe,js_pe:je_pe,n)/10/rho_0*maskZ(is_pe:ie_pe,js_pe:je_pe,nz)
 220   call border_exchg_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,taux(:,:,n)); 
 221   call setcyclic_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,taux(:,:,n))
 222  enddo
 223 
 224  iret=nf_inq_varid(ncid,'tauy',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 225  iret= nf_get_vara_double(ncid,id,(/1,1,1/),(/nx,ny,12/), dat8 )
 226  where (dat8 <= -1e20 ) dat8 = 0
 227  do n=1,12
 228   tauy(is_pe:ie_pe,js_pe:je_pe,n) = dat8(is_pe:ie_pe,js_pe:je_pe,n)/10/rho_0*maskZ(is_pe:ie_pe,js_pe:je_pe,nz)
 229   call border_exchg_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,tauy(:,:,n)); 
 230   call setcyclic_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,tauy(:,:,n))
 231  enddo
 232 
 233  ! heat flux and salinity restoring
 234  iret=nf_inq_varid(ncid,'sst_clim',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 235  iret= nf_get_vara_double(ncid,id,(/1,1,1/),(/nx,ny,12/), dat8 )
 236  where (dat8 <= -1e20 ) dat8 = 0
 237  t_clim(is_pe:ie_pe,js_pe:je_pe,:) = dat8(is_pe:ie_pe,js_pe:je_pe,:)
 238 
 239  iret=nf_inq_varid(ncid,'sst_rest',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 240  iret= nf_get_vara_double(ncid,id,(/1,1,1/),(/nx,ny,12/), dat8 )
 241  where (dat8 <= -1e20 ) dat8 = 0
 242  t_rest(is_pe:ie_pe,js_pe:je_pe,:) = dat8(is_pe:ie_pe,js_pe:je_pe,:)*41868.
 243 
 244  iret=nf_inq_varid(ncid,'sss_clim',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 245  iret= nf_get_vara_double(ncid,id,(/1,1,1/),(/nx,ny,12/), dat8 )
 246  where (dat8 <= -1e20 ) dat8 = 0
 247  s_clim(is_pe:ie_pe,js_pe:je_pe,:) = dat8(is_pe:ie_pe,js_pe:je_pe,:)*1000+35
 248 
 249  iret=nf_inq_varid(ncid,'sss_rest',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 250  iret= nf_get_vara_double(ncid,id,(/1,1,1/),(/nx,ny,12/), dat8 )
 251  where (dat8 <= -1e20 ) dat8 = 0
 252  s_rest(is_pe:ie_pe,js_pe:je_pe,:) = dat8(is_pe:ie_pe,js_pe:je_pe,:)/100.0
 253 
 254  iret = nf_close (ncid)
 255 
 256  ! Restoring zone
 257 
 258  iret=nf_open('restoring_zone.cdf',NF_NOWRITE,ncid)
 259  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 260 
 261  iret=nf_inq_varid(ncid,'t_star',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 262  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1,1/),(/i_blk,j_blk,nz,12/), t_star(is_pe:ie_pe,js_pe:je_pe,:,:) )
 263 
 264  iret=nf_inq_varid(ncid,'s_star',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 265  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1,1/),(/i_blk,j_blk,nz,12/), s_star(is_pe:ie_pe,js_pe:je_pe,:,:) )
 266 
 267  iret=nf_inq_varid(ncid,'tscl',id); if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 268  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1,1/),(/i_blk,j_blk,nz,1/), tscl(is_pe:ie_pe,js_pe:je_pe,:) )
 269 
 270  iret = nf_close (ncid)
 271 
 272  if (enable_idemix ) then
 273   open(10,file='tidal_energy.bin',access='direct',recl=4*nx*ny,status='old')
 274   read(10,rec=1) dat4
 275   close(10)
 276   do j=js_pe,je_pe
 277    do i=is_pe,ie_pe
 278       k=max(1,kbot(i,j))
 279       forc_iw_bottom(i,j) = dat4(i,j)*maskW(i,j,k)/rho_0
 280    enddo
 281   enddo
 282   open(10,file='wind_energy.bin',access='direct',recl=4*nx*ny,status='old')
 283   read(10,rec=1) dat4
 284   close(10)
 285   do j=js_pe,je_pe
 286    do i=is_pe,ie_pe
 287       forc_iw_surface(i,j) = dat4(i,j)*maskW(i,j,nz)/rho_0*0.2
 288    enddo
 289   enddo
 290  endif
 291 end subroutine set_initial_conditions
 292 
 293 
 294 
 295 
 296 
 297 subroutine get_periodic_interval(currentTime,cycleLength,recSpacing,nbrec,trec1,trec2,wght1,wght2)
 298  ! interpolation routine taken from mitgcm
 299  implicit none
 300  real*8, intent(in) :: currentTime,recSpacing,cycleLength
 301  integer, intent(in) :: nbrec
 302  real*8, intent(out) :: wght1,wght2
 303  integer, intent(out) :: tRec1,tRec2
 304  real*8 :: locTime,tmpTime
 305  locTime = currentTime - recSpacing*0.5 + cycleLength*( 2 - NINT(currentTime/cycleLength) )
 306  tmpTime = MOD( locTime, cycleLength )
 307  tRec1 = 1 + INT( tmpTime/recSpacing )
 308  tRec2 = 1 + MOD( tRec1, nbRec )
 309  wght2 = ( tmpTime - recSpacing*(tRec1 - 1) )/recSpacing
 310  wght1 = 1d0 - wght2
 311 end subroutine
 312 
 313 
 314 
 315 
 316 
 317 subroutine set_forcing
 318  use main_module   
 319  use config_module   
 320  use eke_module   
 321  use tke_module   
 322  use idemix_module   
 323  use isoneutral_module   
 324  implicit none
 325  integer :: i,j,n1,n2
 326  real*8 :: t1,t2,f1,f2,fxa,cp_0 = 3991.86795711963  ! J/kg /K
 327  !   q(Tstar-T),  q     (W/m^2/K)   q T/cp0 /rho0  ( W/m^2  K kg /J  m^3/kg = K m /s  )
 328 
 329  fxa = 365*86400d0
 330  call get_periodic_interval((itt-1)*dt_tracer,fxa,fxa/12.,12,n1,n2,f1,f2)
 331 
 332  do j=js_pe-1,je_pe+1
 333   do i=is_pe-1,ie_pe+1
 334     t1 = (taux(i,j-1,n1)*maskZ(i,j-1,nz)+taux(i,j,n1)*maskZ(i,j,nz)) / (maskZ(i,j,nz)+maskZ(i,j-1,nz)+1d-20)
 335     t2 = (taux(i,j-1,n2)*maskZ(i,j-1,nz)+taux(i,j,n2)*maskZ(i,j,nz)) / (maskZ(i,j,nz)+maskZ(i,j-1,nz)+1d-20)
 336     surface_taux(i,j)=(f1*t1+f2*t2)*maskU(i,j,nz)
 337     t1 = (tauy(i-1,j,n1)*maskZ(i-1,j,nz)+tauy(i,j,n1)*maskZ(i,j,nz)) / (maskZ(i,j,nz)+maskZ(i-1,j,nz)+1d-20)
 338     t2 = (tauy(i-1,j,n2)*maskZ(i-1,j,nz)+tauy(i,j,n2)*maskZ(i,j,nz)) / (maskZ(i,j,nz)+maskZ(i-1,j,nz)+1d-20)
 339     surface_tauy(i,j)=( f1*t1 +f2*t2 )*maskV(i,j,nz)
 340   enddo
 341  enddo
 342 
 343  if (enable_tke ) then
 344   do j=js_pe-onx+1,je_pe+onx
 345    do i=is_pe-onx+1,ie_pe+onx
 346      forc_tke_surface(i,j) = sqrt( (0.5*(surface_taux(i,j)+surface_taux(i-1,j)))**2  &
 347                                   +(0.5*(surface_tauy(i,j)+surface_tauy(i,j-1)))**2 )**(3./2.) 
 348    enddo
 349   enddo
 350  endif
 351 
 352  do j=js_pe,je_pe 
 353   do i=is_pe,ie_pe
 354    forc_temp_surface(i,j)=(f1*t_rest(i,j,n1)+f2*t_rest(i,j,n2))* &
 355                           (f1*t_clim(i,j,n1)+f2*t_clim(i,j,n2)-temp(i,j,nz,tau))*maskT(i,j,nz) /cp_0/rho_0
 356    forc_salt_surface(i,j)=(f1*s_rest(i,j,n1)+f2*s_rest(i,j,n2))* &
 357                           (f1*s_clim(i,j,n1)+f2*s_clim(i,j,n2)-salt(i,j,nz,tau))*maskT(i,j,nz)
 358      if (temp(i,j,nz,tau)*maskT(i,j,nz) <= -1.8 .and. forc_temp_surface(i,j) <=0 ) then ! apply simple ice mask
 359        forc_temp_surface(i,j) = 0.0
 360        forc_salt_surface(i,j) = 0.0
 361      endif
 362   enddo
 363  enddo
 364 
 365  if (enable_tempsalt_sources) then
 366   do j=js_pe,je_pe 
 367    do i=is_pe,ie_pe
 368       temp_source(i,j,:) = maskT(i,j,:)*tscl(i,j,:)*(f1*t_star(i,j,:,n1)+f2*t_star(i,j,:,n2) - temp(i,j,:,tau) )
 369       salt_source(i,j,:) = maskT(i,j,:)*tscl(i,j,:)*(f1*s_star(i,j,:,n1)+f2*s_star(i,j,:,n2) - salt(i,j,:,tau) )
 370    enddo
 371   enddo
 372  endif
 373 
 374 end subroutine set_forcing
 375 
 376 
 377 
 378 
 379 subroutine set_diagnostics
 380  use main_module   
 381  use eke_module   
 382  use tke_module   
 383  use idemix_module   
 384  use isoneutral_module   
 385  implicit none
 386  call register_average('taux','Zonal wind stress','m^2/s','UT',surface_taux,0D0,.false.)
 387  call register_average('tauy','Meridional wind stress','m^2/s','TU',surface_tauy,0D0,.false.)
 388  call register_average('forc_temp_surface','Surface temperature flux','m K/s','TT',forc_temp_surface,0D0,.false.)
 389  call register_average('forc_salt_surface','Surface salinity flux','m g/s kg','TT',forc_salt_surface,0D0,.false.)
 390  if (enable_streamfunction) then
 391    call register_average('psi','Barotropic streamfunction','m^2/s','UU',psi(:,:,tau),0D0,.false.)
 392  else
 393    call register_average('psi','Surface pressure','m^2/s','TT',psi(:,:,tau),0D0,.false.)
 394  endif
 395  call register_average('temp','Temperature','deg C','TTT',0d0,temp(:,:,:,tau),.true.)
 396  call register_average('salt','Salinity','g/kg','TTT',0d0,salt(:,:,:,tau),.true.)
 397  call register_average('u','Zonal velocity','m/s','UTT',0d0,u(:,:,:,tau),.true.)
 398  call register_average('v','Meridional velocity','m/s','TUT',0d0,v(:,:,:,tau),.true.)
 399  call register_average('w','Vertical velocity','m/s','TTU',0d0,w(:,:,:,tau),.true.)
 400  call register_average('kappaH','Vertical diffusivity','m^2/s','TTU',0d0,kappaH,.true.)
 401  if (enable_skew_diffusion)  then
 402    call register_average('B1_gm','Zonal component of GM streamfct.','m^2/s','TUT',0d0,B1_gm,.true.)
 403    call register_average('B2_gm','Meridional component of GM streamfct.','m^2/s','UTT',0d0,B2_gm,.true.)
 404  endif
 405  if (enable_tke)  then
 406    call register_average('TKE','Turbulent kinetic energy','m^2/s^2','TTU',0d0,tke(:,:,:,tau),.true.)
 407    call register_average('Prandtl','Prandtl number',' ','TTU',0d0,Prandtlnumber,.true.)
 408    call register_average('mxl','Mixing length',' ','TTU',0d0,mxl,.true.)
 409    call register_average('tke_diss','Dissipation of TKE','m^2/s^3','TTU',0d0,tke_diss,.true.)
 410    call register_average('forc_tke_surface','TKE surface forcing','m^3/s^2','TT',forc_tke_surface,0D0,.false.)
 411    call register_average('tke_surface_corr','TKE surface flux correction','m^3/s^2','TT',tke_surf_corr,0D0,.false.)
 412  endif
 413  if (enable_idemix)  then
 414    call register_average('E_iw','Internal wave energy','m^2/s^2','TTU',0d0,e_iw(:,:,:,tau),.true.)
 415    call register_average('forc_iw_surface','IW surface forcing','m^3/s^2','TT',forc_iw_surface,0D0,.false.)
 416    call register_average('forc_iw_bottom','IW bottom forcing','m^3/s^2','TT',forc_iw_bottom,0D0,.false.)
 417    call register_average('iw_diss','Dissipation of E_iw','m^2/s^3','TTU',0d0,iw_diss,.true.)
 418  endif
 419  if (enable_eke)  then
 420   call register_average('EKE','Eddy energy','m^2/s^2','TTU',0d0,eke(:,:,:,tau),.true.)
 421   call register_average('K_gm','Lateral diffusivity','m^2/s','TTU',0d0,K_gm,.true.)
 422   call register_average('eke_diss_tke','Eddy energy dissipation','m^2/s^3','TTU',0d0,eke_diss_tke,.true.)
 423   call register_average('eke_diss_iw','Eddy energy dissipation','m^2/s^3','TTU',0d0,eke_diss_iw,.true.)
 424  endif
 425 
 426 end subroutine set_diagnostics
 427 
 428 subroutine set_particles
 429 end subroutine set_particles