1 !=======================================================================
   2 ! global 4x4 deg model with 15 levels
   3 !======================================================================= 
   4 
   5 
   6 module config_module
   7  ! use this module only locally in this file
   8  implicit none
   9  real*8, allocatable :: t_star(:,:,:)
  10  real*8, allocatable :: s_star(:,:,:)
  11  real*8, allocatable :: qnec(:,:,:)
  12  real*8, allocatable :: qnet(:,:,:)
  13  real*8, allocatable :: taux(:,:,:)
  14  real*8, allocatable :: tauy(:,:,:)
  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 eke_module   
  25  use tke_module   
  26  use idemix_module   
  27  use isoneutral_module   
  28  use diagnostics_module   
  29  implicit none
  30   nx   = 90
  31   ny   = 40
  32   nz   = 15
  33   dt_mom    = 1800.0
  34   dt_tracer = 86400.!dt_mom
  35 
  36   coord_degree     = .true.
  37   enable_cyclic_x  = .true.
  38 
  39   runlen = 365*86400*100.
  40 !  runlen = 365*86400/12.0
  41   enable_diag_ts_monitor = .true.; ts_monint = 365*86400./24.
  42   enable_diag_snapshots  = .true.; snapint  =  365*86400. /24.0
  43   enable_diag_overturning= .true.; overint  =  365*86400./24.0; overfreq = dt_tracer
  44   enable_diag_energy     = .true.; 
  45   energint = 365*86400./24.
  46   energfreq = 86400
  47   enable_diag_averages   = .true.
  48   aveint  = 365*86400.*10
  49   avefreq = 86400
  50 
  51 
  52   congr_epsilon = 1e-8
  53   congr_max_iterations = 20000
  54   enable_streamfunction = .true.
  55 !  enable_congrad_verbose = .true.
  56 
  57   enable_neutral_diffusion = .true.; 
  58       K_iso_0 = 1000.0
  59       K_iso_steep = 1000.0
  60       iso_dslope=4./1000.0
  61       iso_slopec=1./1000.0
  62   enable_skew_diffusion = .true.
  63 
  64   enable_hor_friction  = .true.; A_h = (4*degtom)**3*2e-11
  65   enable_hor_friction_cos_scaling = .true.; hor_friction_cosPower=1
  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 = .true.
  75 
  76   enable_eke = .true.
  77   eke_k_max  = 1e4
  78   eke_c_k    = 0.4
  79   eke_c_eps  = 0.5
  80   eke_cross  = 2.
  81   eke_crhin  = 1.0
  82   eke_lmin   = 100.0
  83   enable_eke_superbee_advection = .true.
  84 
  85   enable_idemix = .true.
  86   enable_idemix_hor_diffusion = .true.;
  87   enable_eke_diss_surfbot = .true.
  88   eke_diss_surfbot_frac = 0.2 ! fraction which goes into bottom
  89   enable_idemix_superbee_advection = .true.
  90 
  91   eq_of_state_type = 5
  92 end subroutine set_parameter
  93 
  94 
  95 subroutine set_grid
  96   use main_module   
  97   use config_module   
  98   implicit none
  99   real*8 :: ddz(nz)
 100   ddz = (/50.,70.,100.,140.,190.,240.,290.,340.,390.,440.,490.,540.,590.,640.,690./)
 101   dzt=ddz(nz:1:-1)
 102   dxt = 4.0
 103   dyt = 4.0
 104   y_origin = -76.0
 105   x_origin = 4.0
 106 end subroutine set_grid
 107 
 108 
 109 subroutine set_coriolis
 110  use main_module   
 111  use config_module   
 112  implicit none
 113  integer :: j
 114  do j=js_pe-onx,je_pe+onx
 115    coriolis_t(:,j) = 2*omega*sin( yt(j)/180.*pi ) 
 116  enddo
 117 end subroutine set_coriolis
 118 
 119 
 120 
 121 
 122 subroutine set_topography
 123  use main_module   
 124  implicit none
 125  integer :: i,j,k,kk
 126  real*4 :: lev_salt(nx,ny,nz),bathy(nx,ny)
 127  kbot=0
 128 
 129  open(10,file='bathymetry.bin',access='direct',recl=4*nx*ny,status='old')
 130  read(10,rec=1) bathy
 131  close(10)
 132 
 133  open(10,file='lev_s.bin',access='direct',recl=4*nx*ny,status='old')
 134  kk=nz
 135  do k=1,nz
 136   read(10,rec=kk) lev_salt(:,:,k)
 137   kk=kk-1
 138  enddo
 139  close(10)
 140 
 141  do i=is_pe,ie_pe
 142   do j=js_pe,je_pe
 143    do k=nz,1,-1
 144     if (lev_salt(i,j,k) .ne. 0.0 ) kbot(i,j)=k
 145    enddo
 146    if (bathy(i,j) == 0.0) kbot(i,j) =0
 147   enddo
 148  enddo
 149  where ( kbot == nz) kbot = 0
 150  !where (bathy == 0.0) kbot(1:nx,1:ny) =0
 151 end subroutine set_topography
 152 
 153 
 154 subroutine set_initial_conditions
 155  use main_module   
 156  use config_module   
 157  use eke_module   
 158  use tke_module   
 159  use idemix_module   
 160  implicit none
 161  integer :: i,j,k,kk,n
 162  real*4 :: dat4(nx,ny)
 163  include "netcdf.inc"
 164  integer :: iret, ncid,id
 165  real*8 :: fxa,fxb
 166 
 167  allocate( s_star(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); s_star=0
 168  allocate( t_star(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); t_star=0
 169  allocate( qnec(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ) ; qnec=0
 170  allocate( qnet(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); qnet=0
 171  allocate( taux(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); taux=0
 172  allocate( tauy(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); tauy=0
 173 
 174  iret=nf_open('ecmwf_4deg_monthly.cdf',NF_NOWRITE,ncid) ! ECMWF monthly heat flux forcing
 175  !iret=nf_open('ecmwf_4deg.cdf',NF_NOWRITE,ncid) ! annual mean
 176  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 177  iret=nf_inq_varid(ncid,'Q3',id)
 178  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,12/),qnec(is_pe:ie_pe,js_pe:je_pe,1:12))
 179  iret = nf_close (ncid)
 180  do n=1,12
 181   qnec(:,:,n)=qnec(:,:,n)*maskT(:,:,nz)
 182  enddo
 183 
 184 
 185  ! use NCEP net heat flux instead of ECMWF
 186  open(10,file='ncep_qnet.bin',access='direct',recl=4*nx*ny,status='old')
 187  do n=1,12
 188   read(10,rec=n) dat4
 189   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)
 190  enddo
 191  close(10)
 192 
 193  fxa = 0.; fxb = 0.
 194  do n=1,12
 195   do i=is_pe,ie_pe
 196    do j=js_pe,je_pe
 197     fxa = fxa + qnet(i,j,n)*area_t(i,j)
 198     fxb = fxb + area_t(i,j)
 199    enddo
 200   enddo
 201  enddo
 202 
 203  call global_sum(fxa); call global_sum(fxb)
 204  if (my_pe==0) print*,' removing a heat flux imbalance of ',fxa/fxb,'W/m^2'
 205  do n=1,12
 206   qnet(:,:,n) = (qnet(:,:,n) - fxa/fxb)*maskT(:,:,nz)
 207  enddo
 208 
 209  ! use Trenberth wind stress from MITgcm instead of ECMWF (also contained in ecmwf_4deg.cdf)
 210  open(10,file='trenberth_taux.bin',access='direct',recl=4*nx*ny,status='old')
 211  open(20,file='trenberth_tauy.bin',access='direct',recl=4*nx*ny,status='old')
 212  do n=1,12
 213   read(10,rec=n) dat4
 214   taux(is_pe:ie_pe,js_pe:je_pe,n) = dat4(is_pe:ie_pe,js_pe:je_pe)/rho_0
 215   read(20,rec=n) dat4
 216   tauy(is_pe:ie_pe,js_pe:je_pe,n) = dat4(is_pe:ie_pe,js_pe:je_pe)/rho_0
 217   call border_exchg_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,taux(:,:,n)) 
 218   call setcyclic_xy   (is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,taux(:,:,n))
 219   call border_exchg_xy(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,tauy(:,:,n)) 
 220   call setcyclic_xy   (is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,tauy(:,:,n))
 221  enddo
 222  close(10); close(20)
 223  ! check for special values
 224  where( taux < -99.9) taux = 0
 225  where( tauy < -99.9) tauy = 0
 226 
 227  ! initial conditions for T and S
 228  open(10,file='lev_t.bin',access='direct',recl=4*nx*ny,status='old')
 229  open(20,file='lev_s.bin',access='direct',recl=4*nx*ny,status='old')
 230  kk=nz
 231  do k=1,nz
 232   read(10,rec=kk) dat4
 233   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)
 234   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)
 235   read(20,rec=kk) dat4
 236   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)
 237   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)
 238   kk=kk-1
 239  enddo
 240  close(10); close(20)
 241 
 242  ! SST and SSS
 243  open(10,file='lev_sst.bin',access='direct',recl=4*nx*ny,status='old')
 244  open(20,file='lev_sss.bin',access='direct',recl=4*nx*ny,status='old')
 245  do n=1,12
 246   read(10,rec=n) dat4
 247   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)
 248   read(20,rec=n) dat4
 249   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)
 250  enddo
 251  close(10); close(20)
 252 
 253 
 254  if (enable_idemix ) then
 255   open(10,file='tidal_energy.bin',access='direct',recl=4*nx*ny,status='old')
 256   read(10,rec=1) dat4
 257   close(10)
 258   do j=js_pe,je_pe
 259    do i=is_pe,ie_pe
 260       k=max(1,kbot(i,j))
 261       forc_iw_bottom(i,j) = dat4(i,j)*maskW(i,j,k)/rho_0
 262    enddo
 263   enddo
 264   open(10,file='wind_energy.bin',access='direct',recl=4*nx*ny,status='old')
 265   read(10,rec=1) dat4
 266   close(10)
 267   do j=js_pe,je_pe
 268    do i=is_pe,ie_pe
 269       forc_iw_surface(i,j) = dat4(i,j)*maskW(i,j,nz)/rho_0*0.2
 270    enddo
 271   enddo
 272  endif
 273 end subroutine set_initial_conditions
 274 
 275 
 276 
 277 
 278 
 279 subroutine get_periodic_interval(currentTime,cycleLength,recSpacing,nbrec,trec1,trec2,wght1,wght2)
 280  ! interpolation routine taken from mitgcm
 281  implicit none
 282  real*8, intent(in) :: currentTime,recSpacing,cycleLength
 283  integer, intent(in) :: nbrec
 284  real*8, intent(out) :: wght1,wght2
 285  integer, intent(out) :: tRec1,tRec2
 286  real*8 :: locTime,tmpTime
 287  locTime = currentTime - recSpacing*0.5 + cycleLength*( 2 - NINT(currentTime/cycleLength) )
 288  tmpTime = MOD( locTime, cycleLength )
 289  tRec1 = 1 + INT( tmpTime/recSpacing )
 290  tRec2 = 1 + MOD( tRec1, nbRec )
 291  wght2 = ( tmpTime - recSpacing*(tRec1 - 1) )/recSpacing
 292  wght1 = 1d0 - wght2
 293 end subroutine
 294 
 295 
 296 subroutine set_forcing
 297  use main_module   
 298  use config_module   
 299  use eke_module   
 300  use tke_module   
 301  use idemix_module   
 302  implicit none
 303  integer :: i,j,n1,n2
 304  real*8 :: t_rest= 30*86400, cp_0 = 3991.86795711963  ! J/kg /K
 305  real*8 :: fxa,qqnet,qqnec,f1,f2
 306 
 307  fxa = 365*86400d0
 308  call get_periodic_interval((itt-1)*dt_tracer,fxa,fxa/12.,12,n1,n2,f1,f2)
 309 
 310  ! linearly interpolate wind stress and shift from MITgcm U/V grid to this grid
 311  do j=js_pe-onx,je_pe+onx-1
 312   do i=is_pe-onx,ie_pe+onx-1
 313    surface_taux(i,j) = f1*taux(i+1,j,n1) + f2*taux(i+1,j,n2)
 314    surface_tauy(i,j) = f1*tauy(i,j+1,n1) + f2*tauy(i,j+1,n2)
 315   enddo
 316  enddo
 317 
 318  if (enable_tke ) then
 319   do j=js_pe-onx+1,je_pe+onx
 320    do i=is_pe-onx+1,ie_pe+onx
 321      forc_tke_surface(i,j) = sqrt( (0.5*(surface_taux(i,j)+surface_taux(i-1,j)))**2  &
 322                                   +(0.5*(surface_tauy(i,j)+surface_tauy(i,j-1)))**2 )**(3./2.) 
 323    enddo
 324   enddo
 325  endif
 326 
 327  !  W/m^2 K kg/J m^3/kg = K m/s 
 328  do j=js_pe,je_pe 
 329   do i=is_pe,ie_pe
 330    fxa  =  f1*t_star(i,j,n1)+f2*t_star(i,j,n2)
 331    qqnec =  f1*qnec(i,j,n1)+f2*qnec(i,j,n2)
 332    qqnet =  f1*qnet(i,j,n1)+f2*qnet(i,j,n2)
 333    forc_temp_surface(i,j)=(qqnet+ qqnec*(fxa -temp(i,j,nz,tau)) )*maskT(i,j,nz)/cp_0/rho_0
 334 
 335    fxa  =  f1*s_star(i,j,n1)+f2*s_star(i,j,n2)
 336    forc_salt_surface(i,j)=dzt(nz)/t_rest*(fxa-salt(i,j,nz,tau))*maskT(i,j,nz)
 337   enddo
 338  enddo
 339 
 340  ! apply simple ice mask  
 341  do j=js_pe,je_pe 
 342   do i=is_pe,ie_pe
 343    if (temp(i,j,nz,tau)*maskT(i,j,nz) <= -1.8 .and. forc_temp_surface(i,j) <=0 ) then
 344        forc_temp_surface(i,j) = 0.0
 345        forc_salt_surface(i,j) = 0.0
 346    endif
 347   enddo
 348  enddo
 349 
 350 end subroutine set_forcing
 351 
 352 
 353 
 354 
 355 
 356 
 357 subroutine set_diagnostics
 358  use main_module   
 359  use eke_module   
 360  use tke_module   
 361  use idemix_module   
 362  use isoneutral_module   
 363  implicit none
 364  call register_average('taux','Zonal wind stress','m^2/s','UT',surface_taux,0D0,.false.)
 365  call register_average('tauy','Meridional wind stress','m^2/s','TU',surface_tauy,0D0,.false.)
 366  call register_average('forc_temp_surface','Surface temperature flux','m K/s','TT',forc_temp_surface,0D0,.false.)
 367  call register_average('forc_salt_surface','Surface salinity flux','m g/s kg','TT',forc_salt_surface,0D0,.false.)
 368  if (enable_streamfunction) then
 369    call register_average('psi','Barotropic streamfunction','m^2/s','UU',psi(:,:,tau),0D0,.false.)
 370  else
 371    call register_average('psi','Surface pressure','m^2/s','TT',psi(:,:,tau),0D0,.false.)
 372  endif
 373  call register_average('temp','Temperature','deg C','TTT',0d0,temp(:,:,:,tau),.true.)
 374  call register_average('salt','Salinity','g/kg','TTT',0d0,salt(:,:,:,tau),.true.)
 375  call register_average('u','Zonal velocity','m/s','UTT',0d0,u(:,:,:,tau),.true.)
 376  call register_average('v','Meridional velocity','m/s','TUT',0d0,v(:,:,:,tau),.true.)
 377  call register_average('w','Vertical velocity','m/s','TTU',0d0,w(:,:,:,tau),.true.)
 378  call register_average('kappaH','Vertical diffusivity','m^2/s','TTU',0d0,kappaH,.true.)
 379  if (enable_skew_diffusion)  then
 380    call register_average('B1_gm','Zonal component of GM streamfct.','m^2/s','TUT',0d0,B1_gm,.true.)
 381    call register_average('B2_gm','Meridional component of GM streamfct.','m^2/s','UTT',0d0,B2_gm,.true.)
 382  endif
 383  if (enable_tke)  then
 384    call register_average('TKE','Turbulent kinetic energy','m^2/s^2','TTU',0d0,tke(:,:,:,tau),.true.)
 385    call register_average('Prandtl','Prandtl number',' ','TTU',0d0,Prandtlnumber,.true.)
 386    call register_average('mxl','Mixing length',' ','TTU',0d0,mxl,.true.)
 387    call register_average('tke_diss','Dissipation of TKE','m^2/s^3','TTU',0d0,tke_diss,.true.)
 388    call register_average('forc_tke_surface','TKE surface forcing','m^3/s^2','TT',forc_tke_surface,0D0,.false.)
 389    call register_average('tke_surface_corr','TKE surface flux correction','m^3/s^2','TT',tke_surf_corr,0D0,.false.)
 390  endif
 391  if (enable_idemix)  then
 392    call register_average('E_iw','Internal wave energy','m^2/s^2','TTU',0d0,e_iw(:,:,:,tau),.true.)
 393    call register_average('forc_iw_surface','IW surface forcing','m^3/s^2','TT',forc_iw_surface,0D0,.false.)
 394    call register_average('forc_iw_bottom','IW bottom forcing','m^3/s^2','TT',forc_iw_bottom,0D0,.false.)
 395    call register_average('iw_diss','Dissipation of E_iw','m^2/s^3','TTU',0d0,iw_diss,.true.)
 396  endif
 397  if (enable_eke)  then
 398   call register_average('EKE','Eddy energy','m^2/s^2','TTU',0d0,eke(:,:,:,tau),.true.)
 399   call register_average('K_gm','Lateral diffusivity','m^2/s','TTU',0d0,K_gm,.true.)
 400   call register_average('eke_diss_tke','Eddy energy dissipation','m^2/s^3','TTU',0d0,eke_diss_tke,.true.)
 401   call register_average('eke_diss_iw','Eddy energy dissipation','m^2/s^3','TTU',0d0,eke_diss_iw,.true.)
 402  endif
 403 
 404 end subroutine set_diagnostics
 405 
 406 
 407 subroutine set_particles
 408 end subroutine set_particles