1 !=======================================================================
   2 ! global 2x2 deg model with 45 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 :: qsol(:,:,:)
  14  real*8, allocatable :: taux(:,:,:),tauy(:,:,:)
  15  real*8, allocatable :: divpen_shortwave(:)
  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  implicit none
  26   nx   = 128
  27   ny   = 64
  28   nz   = 45
  29   dt_mom    = 3600.0
  30   dt_tracer = 3600.0*5
  31 
  32   coord_degree     = .true.
  33   enable_cyclic_x  = .true.
  34 
  35   runlen = 365*86400.*50 
  36   enable_diag_ts_monitor = .true.; ts_monint = 365*86400./24.
  37   enable_diag_snapshots  = .true.; snapint  =  365*86400.*10
  38   enable_diag_energy     = .true.; energint  = 365*86400.; 
  39   energfreq = dt_tracer*10
  40   enable_diag_averages   = .true.; aveint  = 365*86400.*10
  41   avefreq = dt_tracer*10
  42   enable_diag_overturning  = .true.; overint  = 365*86400.
  43   avefreq = dt_tracer*10
  44 
  45   congr_epsilon = 1e-8
  46   congr_max_iterations = 20000
  47   enable_streamfunction = .true.
  48 
  49   !enable_hor_diffusion = .true.;  K_h = 2000.0
  50   enable_neutral_diffusion = .true.; 
  51       K_iso_0 = 1000.0
  52       K_iso_steep = 500.0
  53       iso_dslope=0.001
  54       iso_slopec=0.001
  55   enable_skew_diffusion = .true.
  56 
  57   enable_hor_friction = .true.; A_h = (2.8*degtom)**3*2e-11    
  58   enable_hor_friction_cos_scaling = .true. 
  59   hor_friction_cosPower = 1
  60   enable_tempsalt_sources = .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 = .true.
  70 
  71   K_gm_0 = 1000.0
  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   enable_eke_superbee_advection = .true.
  80   enable_eke_isopycnal_diffusion = .true.
  81 
  82   enable_idemix = .true.
  83   enable_idemix_hor_diffusion = .true.; 
  84   enable_eke_diss_surfbot = .true.
  85   eke_diss_surfbot_frac = 0.2 ! fraction which goes into bottom
  86   enable_idemix_superbee_advection = .true.
  87 
  88   eq_of_state_type = 5 
  89 end subroutine set_parameter
  90 
  91 
  92 subroutine set_grid
  93   use main_module   
  94   use config_module   
  95   implicit none
  96   real*4 :: dz4(nz)
  97 
  98   open(10,file='dz.bin',access='direct',recl=4*nz,status='old')
  99   read(10,rec=1) dz4
 100   close(10);
 101   dzt = dz4(nz:1:-1)
 102 
 103   dxt = 2.8125
 104   dyt = 2.8125
 105   y_origin = -90.0+2.8125
 106   x_origin = 0+2.8125
 107 end subroutine set_grid
 108 
 109 
 110 subroutine set_coriolis
 111  use main_module   
 112  use config_module   
 113  implicit none
 114  integer :: j
 115  do j=js_pe-onx,je_pe+onx
 116    coriolis_t(:,j) = 2*omega*sin( yt(j)/180.*pi ) 
 117  enddo
 118 end subroutine set_coriolis
 119 
 120 
 121 
 122 subroutine set_initial_conditions
 123  use main_module   
 124  use config_module   
 125  implicit none
 126  integer :: i,j,k,n
 127  real*4 :: dat4(nx,ny)
 128  include "netcdf.inc"
 129  integer :: iret, ncid,id
 130  real*8 :: pen(0:nz),swarg1,swarg2
 131  real*8 ::  rpart_shortwave  = 0.58
 132  real*8 ::  efold1_shortwave = 0.35
 133  real*8 ::  efold2_shortwave = 23.0
 134 
 135  allocate( s_star(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); s_star=0
 136  allocate( t_star(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); t_star=0
 137  allocate( qnec(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ) ; qnec=0
 138  allocate( qnet(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); qnet=0
 139  allocate( qsol(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); qsol=0
 140  allocate( divpen_shortwave(nz) ); divpen_shortwave=0.0
 141  allocate( taux(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); taux=0
 142  allocate( tauy(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx,12) ); tauy=0
 143 
 144  iret=nf_open('forcing_2deg.cdf',NF_NOWRITE,ncid)
 145  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 146 
 147  iret=nf_inq_varid(ncid,'DQDT2',id)
 148  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))
 149  iret=nf_inq_varid(ncid,'QNET2',id)
 150  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,12/),qnet(is_pe:ie_pe,js_pe:je_pe,1:12))
 151  iret=nf_inq_varid(ncid,'SWF2',id)
 152  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,12/),qsol(is_pe:ie_pe,js_pe:je_pe,1:12))
 153 
 154  iret=nf_inq_varid(ncid,'TAUX2',id)
 155  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,12/),taux(is_pe:ie_pe,js_pe:je_pe,1:12))
 156  iret=nf_inq_varid(ncid,'TAUY2',id)
 157  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,12/),tauy(is_pe:ie_pe,js_pe:je_pe,1:12))
 158 
 159  iret=nf_inq_varid(ncid,'SST2',id)
 160  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,12/),t_star(is_pe:ie_pe,js_pe:je_pe,1:12))
 161  iret=nf_inq_varid(ncid,'SSS2',id)
 162  iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,12/),s_star(is_pe:ie_pe,js_pe:je_pe,1:12))
 163 
 164  iret=nf_inq_varid(ncid,'TEMP2',id)
 165  do k=1,nz
 166   iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,k,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/),temp(is_pe:ie_pe,js_pe:je_pe,k,tau))
 167  enddo
 168  temp(:,:,:,tau) = temp(:,:,:,tau)*maskT
 169  where( temp(:,:,:,tau) <= -1e33) temp(:,:,:,tau)=0.0
 170  call border_exchg_xyz(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,nz,temp(:,:,:,tau)) 
 171  call setcyclic_xyz   (is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,nz,temp(:,:,:,tau))
 172  temp(:,:,:,taum1) = temp(:,:,:,tau)
 173 
 174  iret=nf_inq_varid(ncid,'SALT2',id)
 175  do k=1,nz
 176   iret= nf_get_vara_double(ncid,id,(/is_pe,js_pe,k,1/), (/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/),salt(is_pe:ie_pe,js_pe:je_pe,k,tau))
 177  enddo
 178  salt(:,:,:,tau) = salt(:,:,:,tau)*maskT
 179  where( salt(:,:,:,tau) <= -1e33) salt(:,:,:,tau)=35.0
 180  call border_exchg_xyz(is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,nz,salt(:,:,:,tau)) 
 181  call setcyclic_xyz   (is_pe-onx,ie_pe+onx,js_pe-onx,je_pe+onx,nz,salt(:,:,:,tau))
 182  salt(:,:,:,taum1) = salt(:,:,:,tau)
 183 
 184  iret = nf_close (ncid)
 185 
 186  do n=1,12
 187   qnec(:,:,n)=qnec(:,:,n)*maskT(:,:,nz)
 188   qnet(:,:,n)=-qnet(:,:,n)*maskT(:,:,nz)
 189   qsol(:,:,n)=-qsol(:,:,n)*maskT(:,:,nz)
 190   taux(:,:,n)=taux(:,:,n)*maskU(:,:,nz)/rho_0
 191   tauy(:,:,n)=tauy(:,:,n)*maskV(:,:,nz)/rho_0
 192   t_star(:,:,n)=t_star(:,:,n)*maskT(:,:,nz)
 193   s_star(:,:,n)=s_star(:,:,n)*maskT(:,:,nz)
 194  enddo
 195 
 196  if (enable_idemix ) then
 197   open(10,file='tidal_energy.bin',access='direct',recl=4*nx*ny,status='old')
 198   read(10,rec=1) dat4
 199   close(10)
 200   do j=js_pe,je_pe
 201    do i=is_pe,ie_pe
 202       k=max(1,kbot(i,j))
 203       forc_iw_bottom(i,j) = dat4(i,j)*maskW(i,j,k)/rho_0
 204    enddo
 205   enddo
 206   open(10,file='wind_energy.bin',access='direct',recl=4*nx*ny,status='old')
 207   read(10,rec=1) dat4
 208   close(10)
 209   do j=js_pe,je_pe
 210    do i=is_pe,ie_pe
 211       forc_iw_surface(i,j) = dat4(i,j)*maskW(i,j,nz)/rho_0*0.2
 212    enddo
 213   enddo
 214  endif
 215 
 216 !     Initialize penetration profile for solar radiation
 217 !     and store divergence in divpen
 218 !     note that pen(nz) is set 0.0 instead of 1.0 to compensate for the
 219 !     shortwave part of the total surface flux 
 220   pen = 0.0
 221   do k=1,nz-1
 222      swarg1 = zw(k)/efold1_shortwave
 223      swarg2 = zw(k)/efold2_shortwave
 224      pen(k) = rpart_shortwave*exp(swarg1)+ (1.0-rpart_shortwave)*exp(swarg2)
 225   enddo
 226   do k=1,nz
 227      divpen_shortwave(k) = (pen(k) - pen(k-1))/dzt(k)
 228   enddo
 229 end subroutine set_initial_conditions
 230 
 231 
 232 
 233 
 234 
 235 subroutine get_periodic_interval(currentTime,cycleLength,recSpacing,nbrec,trec1,trec2,wght1,wght2)
 236  ! interpolation routine taken from mitgcm
 237  implicit none
 238  real*8, intent(in) :: currentTime,recSpacing,cycleLength
 239  integer, intent(in) :: nbrec
 240  real*8, intent(out) :: wght1,wght2
 241  integer, intent(out) :: tRec1,tRec2
 242  real*8 :: locTime,tmpTime
 243  locTime = currentTime - recSpacing*0.5 + cycleLength*( 2 - NINT(currentTime/cycleLength) )
 244  tmpTime = MOD( locTime, cycleLength )
 245  tRec1 = 1 + INT( tmpTime/recSpacing )
 246  tRec2 = 1 + MOD( tRec1, nbRec )
 247  wght2 = ( tmpTime - recSpacing*(tRec1 - 1) )/recSpacing
 248  wght1 = 1d0 - wght2
 249 end subroutine
 250 
 251 
 252 
 253 
 254 
 255 
 256 
 257 subroutine set_forcing
 258  use main_module   
 259  use config_module   
 260  implicit none
 261  integer :: i,j,k,n1,n2
 262  real*8 :: t_rest= 30*86400, cp_0 = 3991.86795711963  ! J/kg /K
 263  real*8 :: f1,f2,fxa,qqnet,qqnec,ice(is_pe-onx:ie_pe+onx,js_pe-onx:je_pe+onx)
 264 
 265  fxa = 365*86400d0
 266  call get_periodic_interval((itt-1)*dt_tracer,fxa,fxa/12.,12,n1,n2,f1,f2)
 267 
 268  ! linearly interpolate wind stress 
 269  do j=js_pe-onx,je_pe+onx-1
 270   do i=is_pe-onx,ie_pe+onx-1
 271    surface_taux(i,j) = f1*taux(i,j,n1) + f2*taux(i,j,n2)
 272    surface_tauy(i,j) = f1*tauy(i,j,n1) + f2*tauy(i,j,n2)
 273   enddo
 274  enddo
 275 
 276  if (enable_tke ) then
 277   do j=js_pe-onx+1,je_pe+onx
 278    do i=is_pe-onx+1,ie_pe+onx
 279      forc_tke_surface(i,j) = sqrt( (0.5*(surface_taux(i,j)+surface_taux(i-1,j)))**2  &
 280                                   +(0.5*(surface_tauy(i,j)+surface_tauy(i,j-1)))**2 )**(3./2.) 
 281    enddo
 282   enddo
 283  endif
 284 
 285  !  W/m^2 K kg/J m^3/kg = K m/s 
 286  do j=js_pe,je_pe 
 287   do i=is_pe,ie_pe
 288    fxa   =  f1*t_star(i,j,n1)+f2*t_star(i,j,n2)
 289    qqnec =  f1*qnec(i,j,n1)+f2*qnec(i,j,n2)
 290    qqnet =  f1*qnet(i,j,n1)+f2*qnet(i,j,n2)
 291    forc_temp_surface(i,j)=(qqnet+qqnec*(fxa -temp(i,j,nz,tau)) )*maskT(i,j,nz)/cp_0/rho_0
 292 
 293    fxa  =  f1*s_star(i,j,n1)+f2*s_star(i,j,n2)
 294    forc_salt_surface(i,j)=dzt(nz)/t_rest*(fxa-salt(i,j,nz,tau))*maskT(i,j,nz)
 295   enddo
 296  enddo
 297 
 298  ! apply simple ice mask  
 299  ice = 1.0
 300  do j=js_pe,je_pe 
 301   do i=is_pe,ie_pe
 302    if (temp(i,j,nz,tau)*maskT(i,j,nz) <= -1.8 .and. forc_temp_surface(i,j) <=0 ) then
 303        forc_temp_surface(i,j) = 0.0
 304        forc_salt_surface(i,j) = 0.0
 305        ice(i,j)               = 0.0
 306    endif
 307   enddo
 308  enddo
 309 
 310  ! solar radiation
 311  do k=1,nz
 312   do j=js_pe,je_pe 
 313    do i=is_pe,ie_pe
 314     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
 315    enddo
 316   enddo
 317  enddo
 318 end subroutine set_forcing
 319 
 320 
 321 
 322 
 323 
 324 
 325 subroutine set_topography
 326  use main_module   
 327  implicit none
 328  include "netcdf.inc"
 329  integer :: iret, ncid,id, i,j,k
 330  real*8 :: bathy(nx,ny)
 331 
 332  kbot=0
 333 
 334  iret=nf_open('topo_2deg.cdf',NF_NOWRITE,ncid)
 335  if (iret /=0 ) then; print*,nf_strerror(iret) ; stop; endif
 336  iret=nf_inq_varid(ncid,'TOPO3',id)
 337  iret= nf_get_vara_double(ncid,id,(/1,1/), (/nx,ny/),bathy)
 338  iret = nf_close (ncid)
 339 
 340  do j=js_pe,je_pe
 341   do i=is_pe,ie_pe
 342    if (bathy(i,j) >= 0.0) then
 343      kbot(i,j) = 0  
 344    else if (bathy(i,j) <= zw(1) ) then
 345     kbot(i,j) = 1  
 346    else
 347     k= minloc( (zw-bathy(i,j))**2,1)-1
 348     kbot(i,j) = max(1,min(nz,k))
 349    endif
 350   enddo
 351  enddo
 352 end subroutine set_topography
 353 
 354 
 355 
 356 
 357 
 358 
 359 
 360 
 361 subroutine set_diagnostics
 362  use main_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('Nsqr','Square of stability frequency','1/s^2','TTU',0d0,Nsqr(:,:,:,tau),.true.)
 379  call register_average('Hd','Dynamic enthalpy','m^2/s^2','TTT',0d0,Hd(:,:,:,tau),.true.)
 380 
 381  call register_average('K_diss_v','Dissipation by vertical friction','m^2/s^3','TTU',0d0,K_diss_v,.true.)
 382  call register_average('K_diss_h','Dissipation by lateral friction','m^2/s^3','TTU',0d0,K_diss_h,.true.)
 383  call register_average('K_diss_bot','Dissipation by bottom friction','m^2/s^3','TTU',0d0,K_diss_bot,.true.)
 384  call register_average('P_diss_v','Dissipation by vertical mixing','m^2/s^3','TTU',0d0,P_diss_v,.true.)
 385  call register_average('P_diss_nonlin','Dissipation by nonlinear vert. mix.','m^2/s^3','TTU',0d0,P_diss_nonlin,.true.)
 386  call register_average('P_diss_iso','Dissipation by Redi mixing tensor','m^2/s^3','TTU',0d0,P_diss_iso,.true.)
 387 
 388  call register_average('kappaH','Vertical diffusivity','m^2/s','TTU',0d0,kappaH,.true.)
 389  if (enable_skew_diffusion)  then
 390    call register_average('B1_gm','Zonal component of GM streamfct.','m^2/s','TUT',0d0,B1_gm,.true.)
 391    call register_average('B2_gm','Meridional component of GM streamfct.','m^2/s','UTT',0d0,B2_gm,.true.)
 392  else
 393    call register_average('kappa_gm','Vertical GM viscosity','m^2/s','TTU',0d0,kappa_gm,.true.)
 394    call register_average('K_diss_gm','Dissipation by GM friction','m^2/s^3','TTU',0d0,K_diss_gm,.true.)
 395  endif
 396 
 397  if (enable_tke)  then
 398    call register_average('TKE','Turbulent kinetic energy','m^2/s^2','TTU',0d0,tke(:,:,:,tau),.true.)
 399    call register_average('Prandtl','Prandtl number',' ','TTU',0d0,Prandtlnumber,.true.)
 400    call register_average('mxl','Mixing length',' ','TTU',0d0,mxl,.true.)
 401    call register_average('tke_diss','Dissipation of TKE','m^2/s^3','TTU',0d0,tke_diss,.true.)
 402    call register_average('forc_tke_surface','TKE surface forcing','m^3/s^2','TT',forc_tke_surface,0D0,.false.)
 403    call register_average('tke_surface_corr','TKE surface flux correction','m^3/s^2','TT',tke_surf_corr,0D0,.false.)
 404  endif
 405  if (enable_idemix)  then
 406    call register_average('E_iw','Internal wave energy','m^2/s^2','TTU',0d0,e_iw(:,:,:,tau),.true.)
 407    call register_average('forc_iw_surface','IW surface forcing','m^3/s^2','TT',forc_iw_surface,0D0,.false.)
 408    call register_average('forc_iw_bottom','IW bottom forcing','m^3/s^2','TT',forc_iw_bottom,0D0,.false.)
 409    call register_average('iw_diss','Dissipation of E_iw','m^2/s^3','TTU',0d0,iw_diss,.true.)
 410    call register_average('c0','Vertical IW group velocity','m/s','TTU',0d0,c0,.true.)
 411    call register_average('v0','Horizontal IW group velocity','m/s','TTU',0d0,v0,.true.)
 412  endif
 413  if (enable_eke)  then
 414   call register_average('EKE','Eddy energy','m^2/s^2','TTU',0d0,eke(:,:,:,tau),.true.)
 415   call register_average('K_gm','Lateral diffusivity','m^2/s','TTU',0d0,K_gm,.true.)
 416   call register_average('eke_diss','Eddy energy dissipation','m^2/s^3','TTU',0d0,eke_diss,.true.)
 417   call register_average('L_Rossby','Rossby radius','m','TT',L_rossby,0d0,.false.)
 418   call register_average('L_Rhines','Rhines scale','m','TTU',0d0,L_Rhines,.true.)
 419 
 420  endif
 421 
 422 end subroutine set_diagnostics