mod_main.f90

Go to the documentation of this file.

!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!==============================================================================|
!   GLOBAL LIMITS AND ARRAY SIZING PARAMETERS                                  !
!==============================================================================|

MODULE lims
   USE mod_prec  
   IMPLICIT NONE
   SAVE
 
   INTEGER ngl                !!GLOBAL NUMBER OF ELEMENTS
   INTEGER mgl                !!GLOBAL NUMBER OF NODES
   INTEGER numqbc_gl          !!GLOBAL NUMBER OF FRESHWATER INFLOW NODES (RIVERS)
   INTEGER nobclsf_gl         !!GLOBAL NUMBER OF LONGSHORE FLOW ADJUSTED OPEN BOUNDARY NODES
   INTEGER ndrft_gl           !!GLOBAL NUMBER OF LAGRANGIAN TRACKING PARTICLES 

   INTEGER n                  !!LOCAL NUMBER OF ELEMENTS 
   INTEGER m                  !!LOCAL NUMBER OF NODES
   INTEGER numqbc             !!LOCAL NUMBER OF FRESHWATER INFLOW NODES
   INTEGER nobclsf            !!LOCAL NUMBER OF LONGSHORE FLOW ADJUSTED OPEN BOUNDARY NODES
   INTEGER ndrft              !!LOCAL NUMBER OF LAGRANGIAN TRACKING PARTICLES
   INTEGER nisbce_1           !!LOCAL NUMBER OF ELEMENTS WITH ISBCE = 1
   INTEGER nisbce_2           !!LOCAL NUMBER OF ELEMENTS WITH ISBCE = 2
   INTEGER nisbce_3           !!LOCAL NUMBER OF ELEMENTS WITH ISBCE = 3

   INTEGER kb                 !!NUMBER OF SIGMA LEVELS
   INTEGER kbm1               !!NUMBER OF SIGMA LEVELS-1
   INTEGER kbm2               !!NUMBER OF SIGMA LEVELS-2
   INTEGER myid               !!UNIQUE PROCESSOR ID (1 => NPROCS)
   INTEGER msrid
   INTEGER ioprocid           !!PROCESSOR ID OF THE PRINTER NODE IF PRESENT
!   INTEGER KSL               !!NUMBER OF STANDARD SEA LEVELS 
   INTEGER,POINTER:: nprocs_total !!TOTAL NUMBER OF PROCESSORS
   INTEGER,TARGET :: nprocs       !!NUMBER OF PROCESSORS IN FVCOM
   INTEGER ne                 !!NUMBER OF UNIQUE EDGES (LOCAL DOMAIN ONLY)
   INTEGER ncv                !!NUMBER OF INTERNAL CONTROL VOLUMES (EXTENDED LOCAL ONLY)
   
   INTEGER ncv_i              !!NUMBER OF INTERNAL CONTROL VOLUMES (LOCAL ONLY)
   INTEGER nt                 !!TOTAL OF LOCAL INTERNAL + HALO ELEMENTS
   INTEGER mt                 !!TOTAL OF LOCAL INTERNAL + HALO NODES
   INTEGER mx_nbr_elem        !!MAX NUMBER OF ELEMENTS SURROUNDING A NODE

   REAL(sp) :: memtot,memcnt


END MODULE lims



!==============================================================================|
!   CONTROL VARIABLES                                                          |
!==============================================================================|

MODULE control
   USE mod_time
   USE mod_prec   
   USE netcdf
   IMPLICIT NONE
   SAVE

   ! Stuff set inside FVCOM at compile time:
   LOGICAL :: serial         !!TRUE IF SINGLE PROCESSOR
   LOGICAL ::    msr         !!TRUE IF MASTER PROCESSOR (MYID==1)
   LOGICAL ::    par         !!TRUE IF MULTIPROCESSOR RUN
   LOGICAL ::   ioproc         !!TRUE IF PROCESSOR IS THE IOPROC
   LOGICAL :: in_mpi_io_loop !!TRUE IF OUTNODE IS ALREADY IN THE LOOP
   CHARACTER(LEN=80) prg_name

   INTEGER mpi_fvcom_group         !! PROCESSORS WORKING ON FVCOM
!   INTEGER MPI_OTHER_GROUP        !! ADD OTHER GROUPS HERE 

   ! FVCOM info strings
   CHARACTER(LEN=80) fvcom_version !!STRING DESCRIBING VERSION
   CHARACTER(LEN=80) institution   !!STRING DESCRIBING FVCOM MOTHER SHIP
   CHARACTER(LEN=80) fvcom_website !!STRING DESCRIBING WEBSITE FOR FVCOM INFO 

   ! Set at command line
   CHARACTER(LEN=80) casename      !!LETTER ACRONYM SPECIFYING CASE
   CHARACTER(LEN=80) infofile      !!INFO  FILE            

   ! FOR LOADING RUNFILE
   LOGICAL blank_namelist
   CHARACTER(LEN=80) namelist_name

    !--Parameters in NameList NML_CASE
    CHARACTER(LEN=80) case_title     !!CASE TITLE                                 
    CHARACTER(Len=80) date_format    !!DATE FORMAT: ie 2007-11-19 => YMD 
    CHARACTER(Len=80) timezone       !!Name of TIME ZONE OR NONE
    CHARACTER(Len=80) start_date     !!DATE TO START MODEL
    CHARACTER(Len=80) end_date       !!DATE TO END MODEL 
    CHARACTER(Len=80) date_reference !!USER DEFINED REFERENCE DATE OR DEFAULT
    
    LOGICAL use_real_world_time
    namelist /nml_case/        &
         & case_title,         &
         & timezone,           &
         & date_format,        &
     & date_reference,     &
         & start_date,         &
         & end_date

    !--Paramters in NameList NML_STARTUP
    CHARACTER(LEN=80) startup_type  !!'hotstart' or 'coldstart' 
    CHARACTER(LEN=80) startup_file  !!NAME OF START FILE
    CHARACTER(LEN=80) startup_ts_type  !!TYPE OF TS START
    CHARACTER(LEN=80) startup_uv_type  !!TYPE OF TS START
    CHARACTER(LEN=80) startup_turb_type  !!TYPE OF TS START
    REAL(sp) :: startup_t_vals(2)
    REAL(sp) :: startup_s_vals(2)
    REAL(sp) :: startup_u_vals
    REAL(sp) :: startup_v_vals
    REAL(sp) :: startup_dmax
        
    LOGICAL cmdln_restart

    CHARACTER(LEN=80), parameter :: startup_type_coldstart = 'coldstart'
    CHARACTER(LEN=80), parameter :: startup_type_hotstart = 'hotstart'
    CHARACTER(LEN=80), parameter :: startup_type_crashrestart = 'crashrestart'
    CHARACTER(LEN=80), parameter :: startup_type_forecast = 'forecast'

    LOGICAL forecast_mode

    CHARACTER(LEN=80), parameter :: startup_type_default = 'default'
    CHARACTER(LEN=80), parameter :: startup_type_constant = 'constant'
    CHARACTER(LEN=80), parameter :: startup_type_linear = 'linear'
    CHARACTER(LEN=80), parameter :: startup_type_observed = 'observed'
    CHARACTER(LEN=80), parameter :: startup_type_setvalues = 'set values'


    namelist /nml_startup/     &
         & startup_type,       & ! coldstart, hotstart, crashrestart
         & startup_file,       &
         & startup_uv_type,    &  ! default, set values
         & startup_turb_type,  &  ! default, set values
         & startup_ts_type,    &  ! constant, linear, observed, set values
         & startup_t_vals,     &
         & startup_s_vals,     &
         & startup_u_vals,     &
         & startup_v_vals,     &
         & startup_dmax



    !--Parameters in NameList NML_IO
    CHARACTER(LEN=80) input_dir         !!MAIN   INPUT DIRECTORY
    CHARACTER(LEN=80) output_dir        !!PARENT OUTPUT DIRECTORY
    INTEGER ireport                !!INTERVAL (IINT) FOR REPORTING OF FLOWFIELD STATISTICS
    LOGICAL visit_all_vars         !!SET THE LEVEL OF COMPLEXITY IN VISIT
    LOGICAL wait_for_visit         !!WAIT FOR VISIT TO CONNECT BEFORE BEGINING INTEGRATION
    LOGICAL use_mpi_io_mode        !!TURN ON THE PRINTER NODE FOR MPI JOBS

    namelist /nml_io/         &
         & input_dir,         &
         & output_dir,        &
         & ireport,           &
         & visit_all_vars,    &
         & wait_for_visit,    &
         & use_mpi_io_mode
        
    !--Parameters in NameList NML_INTEGRATION

! FOR EXPLICIT TIME STEP MODEL
    REAL(dp)  extstep_seconds!!EXTERNAL TIME STEP IN SECONDS;
                             !! THIS VALUE WILL BE TRUNCATED TO MICROSECONDS!!
    INTEGER   isplit         !!NUMBER OF ITERATIONS OF EXTERNAL MODE/INTERNAL STEP

! FOR SEMI-IMPLICIT TIME STEP MODEL
    REAL(dp)  intstep_seconds!!INTERNAL TIME STEP IN SECONDS FOR SEMI IMPLICIT;

    INTEGER   iramp          !!RAMP FACTOR USED TO EASE STARTUP = f(IINT)
    REAL(sp)  static_ssh_adj !!WATER LEVEL CONSTANT ADJUSTMENT
    REAL(sp)  min_depth      !!MINIMUM ALLOWABLE DEPTH         

    namelist /nml_integration/  &
         & extstep_seconds,     &
         & isplit,              &
         & iramp,               &
         & min_depth,           &
         & static_ssh_adj
    


    !--Parameters in NameList NML_RESTART
    LOGICAL rst_on                  !!TRUE IF OUTPUT RESART FILES
    CHARACTER(LEN=80) rst_first_out !!DATE OF FIRST RESTART FILE OUTPUT
    CHARACTER(LEN=80) rst_out_interval       !!INTERVAL IN DECIMAL DAYS FOR CREATING A RESTART FILE
    INTEGER  rst_output_stack       !!NUMBER OF TIMESTEPS PER FILE

    CHARACTER(LEN=80) restart_file_name !!NAME OF RESTART FILE
    ! ONLY VALID TILL THE NAME IS INCRIMENTED IF THE STKCNT IS EXCEEDED

    namelist /nml_restart/   &
         & rst_on,           &
         & rst_first_out,    &
         & rst_out_interval, &
         & rst_output_stack

        
    !--Parameters in NameList NML_NETCDF
    LOGICAL nc_on
    CHARACTER(LEN=80) nc_first_out        !! DATE TO START NETCDF OUTPUT
    CHARACTER(LEN=80) nc_out_interval     !! OUTPUT INTERVAL IN DECIMAL DAYS
    INTEGER  nc_output_stack              !! NUMBER OF TIMESTEPS PER FILE
    CHARACTER(LEN=120) nc_subdomain_files !! DOMAIN OUTPUT SPECS
    ! GRID STUFF
    LOGICAL nc_grid_metrics
    LOGICAL nc_file_date
    ! MODEL DATA
    LOGICAL nc_velocity
    LOGICAL nc_salt_temp
    LOGICAL nc_turbulence
    LOGICAL nc_vertical_vel
    LOGICAL nc_average_vel
    LOGICAL nc_vorticity
    LOGICAL nc_nh_qp
    LOGICAL nc_nh_rhs
    LOGICAL nc_ice
    
    ! FORCING DATA
    LOGICAL nc_wind_vel
    LOGICAL nc_wind_stress
    LOGICAL nc_wave_para     
    LOGICAL nc_wave_stress   
    LOGICAL nc_evap_precip
    LOGICAL nc_surface_heat
    LOGICAL nc_groundwater
    LOGICAL nc_bio
    LOGICAL nc_wqm

    CHARACTER(LEN=80) nc_file_name !!NAME OF NC FILE
    ! ONLY VALID TILL THE NAME IS INCRIMENTED IF THE STKCNT IS EXCEEDED

    namelist /nml_netcdf/    &
         & nc_on,            &
         & nc_first_out,     &
         & nc_out_interval,  &
         & nc_output_stack,  &
         & nc_subdomain_files,&
         & nc_grid_metrics,  &
         & nc_file_date,     &
         & nc_velocity,      &
         & nc_salt_temp,     &
         & nc_turbulence,    &
         & nc_average_vel,   &
         & nc_vertical_vel,  &
         & nc_wind_vel,      &
         & nc_wind_stress,   &
         & nc_evap_precip,   &
         & nc_surface_heat,  &
         & nc_groundwater,   &
         & nc_bio,           &
         & nc_wqm,           &
     & nc_vorticity

    !--Parameters in NameList NML_NETCDF_AV
    LOGICAL ncav_on                         !! TURN ON NETCDF AVERAGING
    CHARACTER(LEN=80) ncav_first_out        !! DATE TO START NETCDF AVERAGE OUTPUT
    CHARACTER(LEN=80) ncav_out_interval     !! OUTPUT INTERVAL IN DECIMAL DAYS
    INTEGER  ncav_output_stack              !! NUMBER OF TIMESTEPS PER FILE
    CHARACTER(LEN=120) ncav_subdomain_files !! DOMAIN OUTPUT SPECS
    ! GRID STUFF
    LOGICAL ncav_grid_metrics
    LOGICAL ncav_file_date
    ! MODEL DATA
    LOGICAL ncav_velocity
    LOGICAL ncav_salt_temp
    LOGICAL ncav_turbulence
    LOGICAL ncav_vertical_vel
    LOGICAL ncav_average_vel
    LOGICAL ncav_vorticity
    LOGICAL ncav_nh_qp
    LOGICAL ncav_nh_rhs
    LOGICAL ncav_ice

    ! FORCING DATA
    LOGICAL ncav_wind_vel
    LOGICAL ncav_wind_stress
    LOGICAL ncav_wave_para    
    LOGICAL ncav_wave_stress  
    LOGICAL ncav_evap_precip
    LOGICAL ncav_surface_heat
    LOGICAL ncav_groundwater
    LOGICAL ncav_bio
    LOGICAL ncav_wqm

    CHARACTER(LEN=80) ncav_file_name !!NAME OF NCAV FILE
    ! ONLY VALID TILL THE NAME IS INCRIMENTED IF THE STKCNT IS EXCEEDED

    namelist /nml_netcdf_av/  &
         & ncav_on,           &
         & ncav_first_out,    &
         & ncav_out_interval, &
         & ncav_output_stack, &
         & ncav_subdomain_files,&
         & ncav_grid_metrics,   &
         & ncav_file_date,      &
         & ncav_velocity,       &
         & ncav_salt_temp,      &
         & ncav_turbulence,     &
         & ncav_average_vel,    &
         & ncav_vertical_vel,   &
         
         & ncav_wind_vel,       &
         & ncav_wind_stress,    &
         & ncav_evap_precip,    &
         & ncav_surface_heat,   &
         & ncav_groundwater,    &
         & ncav_bio,            &
         & ncav_wqm,            &
     & ncav_vorticity

    !--Parameters in NameList NML_PHYSICS
    CHARACTER(LEN=80) horizontal_mixing_type
    CHARACTER(LEN=80) horizontal_mixing_file
    CHARACTER(LEN=80) horizontal_mixing_kind
    REAL(sp)  horizontal_mixing_coefficient
    REAL(sp)  horizontal_prandtl_number

! REMOVED HORCON, THIS IS NOW STORED IN THE ARRAYS, NN_HVC AND CC_HVC
!    REAL(SP) HORCON ! FVCOM NAME
    REAL(sp) hprnu  ! FVCOM NAME

    CHARACTER(LEN=80) vertical_mixing_type
    REAL(sp) vertical_mixing_coefficient
    REAL(sp) vertical_prandtl_number

    REAL(sp) umol   ! FVCOM NAME
    REAL(sp) vprnu  ! FVCOM NAME

    CHARACTER(LEN=80) bottom_roughness_kind 
    CHARACTER(LEN=80) bottom_roughness_type 
    CHARACTER(LEN=80) bottom_roughness_file 
    REAL(sp)  bottom_roughness_minimum          !!MINIMUM BOTTOM DRAG COEFFICIENT
    REAL(sp)  bottom_roughness_lengthscale      !!BOTTOM FRICTION DEPTH LENGTH SCALE

    REAL(sp) cbcmin ! FVCOM NAME

    CHARACTER(LEN=80),parameter :: br_orig   = 'orig'
    CHARACTER(LEN=80),parameter :: br_gotm   = 'gotm'

    LOGICAL convective_overturning
    LOGICAL scalar_positivity_control
    LOGICAL barotropic
    CHARACTER(LEN=80) sea_water_density_function
    CHARACTER(LEN=80) baroclinic_pressure_gradient
    LOGICAL temperature_active
    LOGICAL salinity_active
!J. Ge
    ! for tracer advection
    LOGICAL backward_advection
    INTEGER backward_step
!J. Ge
    LOGICAL surface_wave_mixing

    ! AT PRESENT YOU CAN NOT COMPILE WITH WET DRY WITHOUT USING IT!
    LOGICAL wetting_drying_on

    LOGICAL  recalculate_rho_mean
    CHARACTER(LEN=80) interval_rho_mean

    CHARACTER(LEN=80),parameter :: sw_dens1 = 'dens1'
    CHARACTER(LEN=80),parameter :: sw_dens2 = 'dens2'
    CHARACTER(LEN=80),parameter :: sw_dens3 = 'dens3'

    LOGICAL adcor_on
    LOGICAL equator_beta_plane
    LOGICAL noflux_bot_condition

    namelist /nml_physics/               &
         & horizontal_mixing_type,       &
         & horizontal_mixing_file,       &
         & horizontal_mixing_kind,       &
         & horizontal_mixing_coefficient,&
         & horizontal_prandtl_number,    & 
         & vertical_mixing_type,         &
         & vertical_mixing_coefficient,  &
         & vertical_prandtl_number,      &
         & bottom_roughness_type,        &
         & bottom_roughness_kind,        &
         & bottom_roughness_file,        &
         & bottom_roughness_lengthscale, &
         & bottom_roughness_minimum,     &
         & convective_overturning,       &
         & scalar_positivity_control,    &
         & barotropic,                   &
         & baroclinic_pressure_gradient, &
         & sea_water_density_function,   &
         & recalculate_rho_mean,         &
         & interval_rho_mean,            &
         & temperature_active,           &
         & salinity_active,              &
         & surface_wave_mixing,          &
!J. Ge
         & backward_advection,           &
         & backward_step,                &
!J. Ge
         & wetting_drying_on,            &
         & adcor_on,                     &
     & equator_beta_plane,           &
         & noflux_bot_condition


    !--Parameters in NameList NML_SURFACE_FORCING
  
    LOGICAL  wind_on
    CHARACTER(LEN=80) wind_type
    CHARACTER(LEN=80) wind_file
    CHARACTER(LEN=80) wind_kind
    REAL(sp) wind_x
    REAL(sp) wind_y
    
    LOGICAL heating_on
    CHARACTER(LEN=80) heating_type
    CHARACTER(LEN=80) heating_file
    CHARACTER(LEN=80) heating_kind
    REAL(sp) heating_longwave_perctage
    REAL(sp) heating_longwave_lengthscale
    REAL(sp) heating_shortwave_lengthscale
    REAL(sp) heating_radiation
    REAL(sp) heating_netflux

    !FVCOM NAMES
    REAL(sp) rheat ! LONG WAVE PERCENTAGE
    REAL(sp) zeta1 ! LONG WAVE LENGTH
    REAL(sp) zeta2 ! SHORT WAVE LENGTH
    

    LOGICAL precipitation_on
    CHARACTER(LEN=80) precipitation_file    
    CHARACTER(LEN=80) precipitation_kind    
    REAL(sp) :: precipitation_evp
    REAL(sp) :: precipitation_prc

    LOGICAL airpressure_on
    CHARACTER(LEN=80) airpressure_file    
    CHARACTER(LEN=80) airpressure_kind
    REAL(sp) :: airpressure_value

    LOGICAL wave_on               
    CHARACTER(LEN=80) wave_file    
    CHARACTER(LEN=80) wave_kind
    REAL(sp) :: wave_height
    REAL(sp) :: wave_length
    REAL(sp) :: wave_direction
    REAL(sp) :: wave_period
    REAL(sp) :: wave_per_bot
    REAL(sp) :: wave_ub_bot

    ! FORCING KINDS
    CHARACTER(LEN=80),parameter:: cnstnt   = "constant"
    CHARACTER(LEN=80),parameter:: sttc     = "static"
    CHARACTER(LEN=80),parameter:: tmdpndnt = "time dependant"
    CHARACTER(LEN=80),parameter:: prdc     = "periodic"
    CHARACTER(LEN=80),parameter:: vrbl     = "variable"

    CHARACTER(LEN=80), parameter:: speed  = "speed"
    CHARACTER(LEN=80), parameter:: stress = "stress"

!J. Ge
    CHARACTER(LEN=80),parameter:: uniform = "uniform"
    CHARACTER(LEN=80),parameter:: non_uniform = "non-uniform"
!J. Ge

    namelist /nml_surface_forcing/          &
         & wind_on,                         &
         & wind_type,                       &
         & wind_file,                       &
         & wind_kind,                       &
         & wind_x,                          &
         & wind_y,                          &
         & heating_on,                      &
         & heating_type,                    &
         & heating_kind,                    &
         & heating_file,                    &
         & heating_longwave_lengthscale,    &
         & heating_longwave_perctage,       &
         & heating_shortwave_lengthscale,   & 
         & heating_radiation,               &
         & heating_netflux,                 & 
         & precipitation_on,                &
         & precipitation_kind,              &
         & precipitation_file,              &
         & precipitation_prc,               &
         & precipitation_evp,               &
     & airpressure_on,                  &
     & airpressure_kind,                &
     & airpressure_file,                &
     & airpressure_value,               &
         & wave_on,                         &   
         & wave_file,                       &
         & wave_kind,                       &
         & wave_height,                     &
         & wave_length,                     &
         & wave_direction,                  &
         & wave_period,                     &
         & wave_per_bot,                    &
         & wave_ub_bot  

    !--Parameters in NameList NML_RIVER_TYPE
    CHARACTER(LEN=80) river_ts_setting ! METHOD TO CALCULATE T&S AT
    ! THE RIVER MOUTH
    CHARACTER(LEN=80) river_inflow_location ! IS THE FLUX
    CHARACTER(LEN=80) river_kind ! PERIODIC OR VARIABLE
    CHARACTER(LEN=80) river_info_file ! RIVERS INFORMATION FILE
    ! SPECIFIED FOR NODES OR EDGES
    INTEGER river_number ! THE NUMBER OF RIVERS

    namelist /nml_river_type/               &
         & river_number,                    &
         & river_kind,                      &
         & river_ts_setting,                &
         & river_info_file,                 &
         & river_inflow_location
    

    !--Parameters in NameList NML_RIVER
    INTEGER, PARAMETER :: river_char_len=60

    INTEGER, PARAMETER :: max_layers = 100

    CHARACTER(LEN=80) river_name
    CHARACTER(LEN=80) river_file
    INTEGER river_grid_location

    ! ASSUMING MAXIMUM NUMBER OF LAYERS!
    REAL(sp) :: river_vertical_distribution(max_layers)

    namelist /nml_river/                    &
         & river_name,                      &
         & river_file,                      &
         & river_grid_location,             &        
         & river_vertical_distribution

    
    ! THIS TYPE PROVIDES INTERMEDIATE STORAGE FOR RIVER INFO UNTILL
    ! THE RIVER FORCING IS SET UP IN MOD_FORCE
    type river
       CHARACTER(LEN=80) name
       CHARACTER(LEN=80) file
       INTEGER location
       REAL(sp) :: distribution(max_layers)
    end type river

    type(river), Allocatable, DIMENSION(:) :: rivers

    !--Parameters in NameList NML_OPEN_BOUNDARY
    LOGICAL obc_on
    CHARACTER(LEN=80) obc_node_list_file
    LOGICAL obc_elevation_forcing_on
    CHARACTER(LEN=80) obc_elevation_file
    INTEGER obc_ts_type
    LOGICAL obc_temp_nudging
    CHARACTER(LEN=80) obc_temp_file
    REAL(sp) obc_temp_nudging_timescale
    LOGICAL obc_salt_nudging
    CHARACTER(LEN=80) obc_salt_file
    REAL(sp) obc_salt_nudging_timescale
    LOGICAL obc_meanflow
    CHARACTER(LEN=80) obc_meanflow_file
    LOGICAL obc_longshore_flow_on
    CHARACTER(LEN=80) obc_longshore_flow_file
    INTEGER obc_tideout_initial,obc_tideout_interval
    LOGICAL obc_depth_control_on

    namelist /nml_open_boundary_control/    &
         & obc_on,                          &
         & obc_node_list_file,              &
         & obc_elevation_forcing_on,        &
         & obc_elevation_file,              &
         & obc_ts_type,                     &
         & obc_temp_nudging,                &
         & obc_temp_file,                   &
         & obc_temp_nudging_timescale,      &
         & obc_salt_nudging,                &
         & obc_salt_file,                   &
         & obc_salt_nudging_timescale,      &
         & obc_meanflow,                    &
         & obc_meanflow_file,               &
         & obc_tideout_initial,             &
         & obc_tideout_interval,            &
         & obc_longshore_flow_on,           &
         & obc_longshore_flow_file,         &
         & obc_depth_control_on


    !--Parameters in NameList GRID_COORDINATES
    CHARACTER(LEN=80) grid_file
    CHARACTER(LEN=80) grid_file_units
    CHARACTER(LEN=200) projection_reference
    CHARACTER(LEN=80) sigma_levels_file
    CHARACTER(LEN=80) depth_file
    CHARACTER(LEN=80) coriolis_file
    CHARACTER(LEN=80) sponge_file

    LOGICAL use_proj
! THESE ARE THE FVCOM NAMES FOR 
     !SIGMA_LEVELS = KB
     !SIGMA_LAYERS = KBM1 = KB - 1
     !SIGMA_LEVELS - 2 = KBM2

    namelist /nml_grid_coordinates/        &
         & grid_file,                      &
         & grid_file_units,                &
         & projection_reference,           &
         & sigma_levels_file,              &
         & depth_file,                     &
         & coriolis_file,                  &
         & sponge_file

    !--Parameters in NameList NML_GROUNDWATER
    LOGICAL groundwater_on
    CHARACTER(LEN=80) groundwater_kind
    CHARACTER(LEN=80) groundwater_file
    REAL(sp) groundwater_flow
    REAL(sp) groundwater_temp
    LOGICAL groundwater_temp_on
    REAL(sp) groundwater_salt
    LOGICAL groundwater_salt_on



    namelist /nml_groundwater/            &
         & groundwater_on,                &
         & groundwater_temp_on,           &
         & groundwater_salt_on,           &
         & groundwater_kind,              &
         & groundwater_file,              &
         & groundwater_flow,              &
         & groundwater_temp,              &
         & groundwater_salt

    !--Parameters in NameList NML_LAG_PART
    LOGICAL lag_particles_on
    CHARACTER(LEN=80) lag_start_file
    CHARACTER(LEN=80) lag_out_file
    CHARACTER(LEN=80) lag_first_out
    CHARACTER(LEN=80) lag_restart_file
    CHARACTER(LEN=80) lag_out_interval 
    CHARACTER(LEN=80) lag_scal_choice

!!$    LOGICAL LAG_TEMPERATURE
!!$    LOGICAL LAG_SALINITY
!!$    LOGICAL LAG_DENSITY
!!$    LOGICAL LAG_EDDY_VISCOSITY
!!$    LOGICAL LAG_DIFFUSIVITY
  
    namelist /nml_lag/                    &
         & lag_particles_on,              &
         & lag_start_file,                &
         & lag_out_file,                  &
         & lag_first_out,                 &
         & lag_restart_file,              &
         & lag_out_interval,              &
         & lag_scal_choice


!!$         & LAG_TEMPERATURE,               &
!!$         & LAG_SALINITY,                  &
!!$         & LAG_DENSITY,                   &
!!$         & LAG_EDDY_VISCOSITY,            &
!!$         & LAG_DIFFUSIVITY

    !--Parameters in NameList NML_ADDITIONAL_MODELS
!    LOGICAL WATER_QUALITY_MODEL
!    CHARACTER(LEN=80) WATER_QUALITY_MODEL_FILE
    LOGICAL data_assimilation
    CHARACTER(LEN=80) data_assimilation_file
    LOGICAL biological_model
! J. Ge for online biology
    CHARACTER(LEN=80) biological_model_file
!--------------------------

    CHARACTER(LEN=80) startup_bio_type  !!TYPE OF BIO START
    LOGICAL sediment_model
    CHARACTER(LEN=80) sediment_model_file
    CHARACTER(LEN=80) sediment_parameter_type
    CHARACTER(LEN=80) sediment_parameter_file
    CHARACTER(LEN=80) bedflag_type 
    CHARACTER(LEN=80) bedflag_file 

    LOGICAL icing_model
    CHARACTER(LEN=80) icing_forcing_file
    CHARACTER(LEN=80) icing_forcing_kind
    REAL(sp) ::icing_air_temp
    REAL(sp) ::icing_wspd
    
    LOGICAL ice_model
    CHARACTER(LEN=80) ice_forcing_file
    CHARACTER(LEN=80) ice_forcing_kind
    REAL(sp) ::ice_sea_level_pressure
    REAL(sp) ::ice_air_temp
    REAL(sp) ::ice_spec_humidity
    REAL(sp) ::ice_cloud_cover
    REAL(sp) ::ice_shortwave
    CHARACTER(LEN=80) ice_longwave_type
    
    LOGICAL high_latitude_wave

    !ADDITIONAL MODEL DATA
    INTEGER            :: n_sed
    INTEGER, PARAMETER :: n_sed_max = 10
    CHARACTER(LEN=20)  :: sed_names(n_sed_max)
    REAL(sp), ALLOCATABLE :: seddis(:,:)


    REAL(sp), ALLOCATABLE :: biodis(:,:)


    ! FVCOM RUN MODE PARAMETERS
    CHARACTER(LEN=80) fvcom_run_mode
    CHARACTER(LEN=80),parameter :: fvcom_pure_sim         = 'pure sim'
!    CHARACTER(LEN=80),parameter :: FVCOM_NUDGE_AVG_SST = 'nudge avg sst'
!    CHARACTER(LEN=80),parameter :: FVCOM_NUDGE_AVG_TSGRD = 'nudge avg ts'
    CHARACTER(LEN=80),parameter :: fvcom_nudge_oi_assim   = 'nudge or OI assim'   
!    CHARACTER(LEN=80),parameter :: FVCOM_OI_ASSIM      = 'OI ASSIM'
!    CHARACTER(LEN=80),parameter :: FVCOM_KALMAN_RRKF   = 'Kalman RRKF'
    CHARACTER(LEN=80),parameter :: fvcom_rrkf_without_ssa = 'RRKF WITHOUT SSH/SST'
    CHARACTER(LEN=80),parameter :: fvcom_rrkf_with_ssa    = 'RRKF WITH SSH/SST'
    CHARACTER(LEN=80),parameter :: fvcom_enkf_without_ssa = 'ENKF WITHOUT SSH/SST'
    CHARACTER(LEN=80),parameter :: fvcom_enkf_with_ssa    = 'ENKF WITH SSH/SST'
    CHARACTER(LEN=80),parameter :: fvcom_kalman_4         = 'Kalman 4'

    namelist /nml_additional_models/      &
!         & WATER_QUALITY_MODEL,           &
!         & WATER_QUALITY_MODEL_FILE,      &
         & data_assimilation,             &
         & data_assimilation_file,        &
         & biological_model,              &
         & startup_bio_type,              &
!--------- J. Ge for biology --------------
         & biological_model_file,         &
!------------------------------------------
         & sediment_model,                &
         & sediment_model_file,           &
         & sediment_parameter_type,       &
         & sediment_parameter_file,       &
         & bedflag_type,                  &
         & bedflag_file,                  &
         & icing_model,                   &
         & icing_forcing_file,            &
         & icing_forcing_kind,            &
         & icing_air_temp,                &
         & icing_wspd,                    &
         & ice_model,                     &
         & ice_forcing_file,              &
         & ice_forcing_kind,              &
         & ice_sea_level_pressure,        &
         & ice_air_temp,                  &
         & ice_spec_humidity,             &
         & ice_shortwave,                 &
         & ice_longwave_type,                 &
         & ice_cloud_cover,               &
         & high_latitude_wave

    !--Parameters in NameList NML_PROBE
    LOGICAL probes_on
    INTEGER probes_number
    CHARACTER(len=80) probes_file
    
    namelist /nml_probes/      &
         & probes_on,          &
         & probes_number,      &
         & probes_file

    !--Parameters in NameList NML_BOUNDSCHK
    !=> bounds checking
    LOGICAL  :: force_archive = .false.
    LOGICAL  :: boundschk_on = .false.
    INTEGER  :: chk_interval
    REAL(sp) :: veloc_mag_max
    REAL(sp) :: zeta_mag_max
    REAL(sp) :: temp_max
    REAL(sp) :: temp_min
    REAL(sp) :: salt_max
    REAL(sp) :: salt_min
  
    namelist /nml_boundschk/      &
           & boundschk_on,        &
           & chk_interval,        &
           & veloc_mag_max,       &
           & zeta_mag_max,        &
           & temp_max,            &
           & temp_min,            &
           & salt_max,            &
           & salt_min
    

!--Time Variables for FVCOM-----------------------------------------!
    TYPE(time) :: inttime
    TYPE(time) :: exttime
    Type(time) :: rktime
    TYPE(time) :: zerotime

    TYPE(time) :: endtime
    Type(time) :: starttime
    Type(time) :: runfile_starttime
    Type(time) :: referencedate

    TYPE(time) :: delt_rho_mean
    TYPE(time) :: recalc_rho_mean

    ! ZERO PHASE TIME FOR SPECTRAL (NON JULIAN) TIDE
    TYPE(time) :: spectime
    
    REAL(sp)   :: dte        !!EXTERNAL TIME STEP (Seconds)
    REAL(sp)   :: dti        !!INTERNAL TIME STEP (Seconds)
    REAL(sp)   :: ramp

    TYPE(time) :: imdte      !!EXTERNAL TIME STEP 
    TYPE(time) :: imdti      !!INTERNAL TIME STEP 
    
    INTEGER(itime) :: iint     !!INTERNAL TIME STEP ITERATION NUMBER (ISTART => IEND)
    INTEGER(itime) :: iext     !!EXTERNAL TIME STEP ITERATION NUMBER (1 => ISPLIT)
    INTEGER(itime) :: istart   !!STARTING INTERNAL TIME STEP ITERATION NUMBER
    INTEGER(itime) :: iend     !!ENDING INTERNAL TIME STEP ITERATION NUMBER
    INTEGER(itime) :: nsteps   !!Number OF INTERAL TIME STEPS IN SIMULATION
    
    ! Time variables for File IO
    
    INTEGER, PARAMETER:: timeprec = 6
    !! THIS IS THE LENGHT THAT LOOKS NICE FOR GIVEN TIME PREC
    INTEGER, PARAMETER:: datestrlen = 20+timeprec 

!    CHARACTER(LEN=80) :: IO_FILE_DATE
!    CHARACTER(LEN=80) :: IO_timestr
!    real(SP)          :: IO_days
!    integer           :: IO_mjd
!    integer           :: IO_msec
!    integer           :: IO_IINT



!--Constants-------------------------------------------------------------------!

    ! SELECT WHETHER YOU WANT TO MAKE TIME MORE EXACT (SP CAUSES
    ! ROUND OF ERROR BUT IS THE TRADITIONAL PRECISION)
   REAL(dp), PARAMETER, DIMENSION(4) :: alpha_rk = (/0.2500_dp,1.0_dp/3.0_dp,0.5000_dp,1.0_dp/)
!   REAL(SP), PARAMETER, DIMENSION(4) :: ALPHA_RK = (/0.2500_DP,1.0_DP/3.0_DP,0.5000_DP,1.0_DP/)


   REAL(dp), PARAMETER :: grav      = 9.81_sp
   REAL(dp), PARAMETER :: pi        = 3.141592653589793238_dp
   REAL(dp), PARAMETER :: pi2       = 2.0_dp * 3.141592653589793238_dp
   REAL(dp), PARAMETER :: zero      = 0.0_dp 
   REAL(dp), PARAMETER :: one_third = 1.0_dp/3.0_dp 
   REAL(dp), PARAMETER :: rearth    = 6371.0e03_dp   !!Earth Radius in Meters
   REAL(dp), PARAMETER :: deg2rad   = pi2/360.0_dp   !!Radians/Degree
   REAL(dp), PARAMETER :: tpi       = deg2rad*rearth !TPI=pi*rearth/180.=3.14159265/180.0*6371.*1000.
   REAL(dp), PARAMETER :: rofvros   = 0.9775171065_dp!!RATIO OF THE DENSITY OF FRESH AND SEA WATER 1000./1023.   
   real(dp), parameter :: slp0      = 101325.0_sp    !! mean sea surface pressure (Pa)


!--Parameter Controlling Vertical Coordinate Distribution------------
   !----------!
   CHARACTER(LEN=80) :: stype
   CHARACTER(LEN=80), PARAMETER :: stype_uniform= "UNIFORM"
   CHARACTER(LEN=80), PARAMETER :: stype_geometric= "GEOMETRIC"
   CHARACTER(LEN=80), PARAMETER :: stype_tanh= "TANH"
   CHARACTER(LEN=80), PARAMETER :: stype_generalized= "GENERALIZED"
   CHARACTER(LEN=80), PARAMETER :: stype_restart= "RESTART"
!--Sigma Level Parameters for case GEOMETRIC OF UNIFORM------------------------!
   REAL(sp) :: p_sigma      !!PARAMETER CONTROLLING SIGMA LEVEL DISTRIBUTION

!--General Vertical Level Parameters for case TANH ----------------------------!
   REAL(sp) :: du2          !!PARAMETER CONTROLLING LEVEL DISTRIBUTION OF SURFACE
   REAL(sp) :: dl2          !!PARAMETER CONTROLLING LEVEL DISTRIBUTION OF BOTTOM

!--General Vertical Level Parameters for case GENERALIZED ---------------------!
   REAL(sp) :: duu          !!THE UPPER BOUNDARY OF PARALLEL COORDINATE
   REAL(sp) :: dll          !!THE LOWER BOUNDARY OF PARALLEL COORDINATE
   REAL(sp) :: hmin1        !!THE MIN DEPTH AT WHICH THE LAYERS ARE CONSTANT
   
   INTEGER  :: ku           !!THE NUMBERS OF LAYERS ABOVE UPPER BOUNDARY
   INTEGER  :: kl           !!THE NUMBVER OF LAYERS BELOW LOWER BOUNDARY
   
   REAL(sp), ALLOCATABLE :: zku(:)       !!THE DEPTHS OF PARALLEL LAYERS ABOVE UPPER BOUNDARY
   REAL(sp), ALLOCATABLE :: zkl(:)       !!THE DEPTHS OF PARALLEL LAYERS BELOW LOWER BOUNDARY

   REAL(sp) :: hmax         !!GLOBAL MAXIMUM DEPTH IN DEPTH FILE
   REAL(sp) :: hmin         !!GLOBAL MINIMUM DEPTH IN DEPTH FILE


   ! All file identifiers go here
   INTEGER :: ipt                     !! IUNIT FOR LOG FILE OUTPUT
   INTEGER, PUBLIC, PARAMETER :: ipt_base= 7000 ! FOR PAR LOG FILES

   INTEGER, PARAMETER :: testunit = 200 !! TEST TO SEE IF OUTPUT DIR EXISTS/WRITABLE
   INTEGER, PARAMETER :: nmlunit = 10 !! NAMELIST RUN FILE
   INTEGER, PARAMETER :: itsunit=11
   INTEGER, PARAMETER :: obcunit=12
   INTEGER, PARAMETER :: gridunit=13
   INTEGER, PARAMETER :: sigmaunit=14
   INTEGER, PARAMETER :: depthunit=15
   INTEGER, PARAMETER :: coriolisunit=16
   INTEGER, PARAMETER :: spongeunit=17
   INTEGER, PARAMETER :: lsfunit=18
   INTEGER, PARAMETER :: assimunit=19
   INTEGER, PARAMETER :: oiassimunit=23
   INTEGER, PARAMETER :: probeunit=20
   INTEGER, PARAMETER :: julobcunit=21
   INTEGER, PARAMETER :: kfunit=22
   INTEGER, PARAMETER :: nestunit=30
   INTEGER, PARAMETER :: subdunit=31

   INTEGER            :: rivernmlunit


END MODULE control

!==============================================================================|
!==============================================================================|
!==============================================================================|
!     ALL VARS                                                                 |
!     CONATINS:
!      FVCOM VARIABLES
!      N2E3D: simple average from nodes to elements for 3D variables
!      N2E2D: simple average from nodes to elements for 2D variables
!      E2N3D: simple average from elements to nodes for 3D variables
!      E2N2D: simple average from elements to nodes for 2D variables
!==============================================================================|
MODULE all_vars
   USE mod_prec
   USE lims
   USE control 
   IMPLICIT NONE
   SAVE


!--------------------------Temporary Array------------------------------------------!

  INTEGER, ALLOCATABLE :: nvg(:,:)

!--------------------------Global Grid Variables------------------------------------!

!!  REAL(SP), POINTER     :: XG(:)               !!GLOBAL X-COORD AT NODE 
!!  REAL(SP), POINTER     :: YG(:)               !!GLOBAL X-COORD AT NODE 
  REAL(sp), ALLOCATABLE     :: xg(:)               !!GLOBAL X-COORD AT NODE 
  REAL(sp), ALLOCATABLE     :: yg(:)               !!GLOBAL X-COORD AT NODE 
  REAL(sp), ALLOCATABLE :: hg(:)               !!GLOBAL DEPTH AT NODE
!  REAL(SP), ALLOCATABLE :: CORG(:)             !!GLOBAL COORIOLIS AT NODE 
  REAL(sp), ALLOCATABLE :: xcg(:)              !!GLOBAL X-COORD AT FACE CENTER 
  REAL(sp), ALLOCATABLE :: ycg(:)              !!GLOBAL X-COORD AT FACE CENTER 

!--------------------------Grid Metrics---------------------------------------------!

  ! JUST DON'T ALLOCATE THINGS THAT YOU DON'T NEED IN A PARTICULAR
  ! MODULE IT IS A PAIN IN THE NECK TO NO HAVE THE VARIABLE DECLARED
  

!# if !defined(SPHERICAL)
   REAL(sp)              :: vxmin,vymin,vxmax,vymax
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: xm(:)               !!X-COORD AT NODE IN METERS
   REAL(sp), ALLOCATABLE,TARGET :: ym(:)               !!Y-COORD AT NODE IN METERS
   REAL(sp), ALLOCATABLE,TARGET :: xmc(:)              !!X-COORD AT CELL CENTER IN METERS
   REAL(sp), ALLOCATABLE,TARGET :: ymc(:)              !!Y-COORD AT CELL CENTER IN METERS
   REAL(sp), ALLOCATABLE,TARGET :: lon(:)              !!LONGITUDE AT THE NODE
   REAL(sp), ALLOCATABLE,TARGET :: lat(:)              !!LATITUDE AT THE NODE
   REAL(sp), ALLOCATABLE,TARGET :: lonc(:)              !!LONGITUDE AT THE NODE
   REAL(sp), ALLOCATABLE,TARGET :: latc(:)              !!LATITUDE AT THE NODE
   ! VX,VY and XC,YC are used for either meters or spherical depending on
   ! make file option 'SPHERICAL'
   REAL(sp), ALLOCATABLE,TARGET :: vx(:)               !!X-COORD AT GRID POINT
   REAL(sp), ALLOCATABLE,TARGET :: vy(:)               !!Y-COORD AT GRID POINT
   REAL(sp), ALLOCATABLE,TARGET :: xc(:)               !!X-COORD AT FACE CENTER 
   REAL(sp), ALLOCATABLE,TARGET :: yc(:)               !!Y-COORD AT FACE CENTER

   ! NOTES: SHOULD MAKE AN ARRAY TO STORE 1/ART, 1/ART2 and 1/ART2
   ! IT is faster and safer

   REAL(sp), ALLOCATABLE,TARGET :: art(:)              !!AREA OF ELEMENT
   REAL(sp), ALLOCATABLE,TARGET :: art1(:)             !!AREA OF NODE-BASE CONTROl VOLUME
   REAL(sp), ALLOCATABLE,TARGET :: art2(:)             !!AREA OF ELEMENTS AROUND NODE
!--Gravity (vary with latitute) -----------------------------------------------!
   REAL(sp),ALLOCATABLE,TARGET :: grav_n(:),grav_e(:) ! CALCULATED AS A
   ! FUNCTION OF LATITUDE IN SPHERICAL COORDINATES MODEL
   
!----------------Node, Boundary Condition, and Control Volume-----------------------!

   INTEGER, ALLOCATABLE,TARGET :: nv(:,:)             !!NODE NUMBERING FOR ELEMENTS
!   INTEGER, ALLOCATABLE,TARGET :: NVGL(:,:)             !!NODE GLOBAL NUMBERING OF LOCAL ELEMENTS
   INTEGER, ALLOCATABLE,TARGET :: nbe(:,:)            !!INDICES OF ELMNT NEIGHBORS
!   INTEGER, POINTER            :: NBEGL(:,:)          !!GLOBAL INDICES OF LOCAL ELEMENT NEIGHBORS
   INTEGER, ALLOCATABLE,TARGET :: ntve(:)             !! NUMBER OF ELEMENTS SURROUNDING EACH NODE
   INTEGER, ALLOCATABLE,TARGET :: ntsn(:)             !! NUMBER OF NODES SURROUNDING EACH NODE
   INTEGER, ALLOCATABLE,TARGET :: isonb(:)            !!NODE MARKER = 0,1,2
   INTEGER, ALLOCATABLE,TARGET :: isonb_w(:)          !!NODE MARKER = 0,1,2
   INTEGER, ALLOCATABLE,TARGET :: isbc(:)     
   INTEGER, ALLOCATABLE,TARGET :: isbce(:)     
   INTEGER, ALLOCATABLE,TARGET :: iec(:,:)
   INTEGER, ALLOCATABLE,TARGET :: ienode(:,:)
   INTEGER, ALLOCATABLE,TARGET :: nbsn(:,:)            !! INDICES OF NODES SURROUNDING EACH NODE
!   INTEGER, POINTER            :: NBSNGL(:,:)          !! GLOBAL INDICIES OF NODES SURROUNDING EACH LOCAL NODE
   INTEGER, ALLOCATABLE,TARGET :: niec(:,:)
   INTEGER, ALLOCATABLE,TARGET :: ntrg(:)
   INTEGER, ALLOCATABLE,TARGET :: nbve(:,:)            !! INDICIES OF ELEMENTS SURROUNDING EACH NODE
!   INTEGER, POINTER            :: NBVEGL(:,:)          !! GLOBAL INDICIES OF ELEMENTS SURROUNDING EACH LOCAL NODE
   INTEGER, ALLOCATABLE,TARGET :: nbvt(:,:)
   INTEGER, ALLOCATABLE,TARGET :: lisbce_1(:)          !!LIST OF ELEMENTS WITH ISBCE=1
   INTEGER, ALLOCATABLE,TARGET :: lisbce_2(:)          !!LIST OF ELEMENTS WITH ISBCE=2
   INTEGER, ALLOCATABLE,TARGET :: lisbce_3(:)          !!LIST OF ELEMENTS WITH ISBCE=3
   REAL(sp),ALLOCATABLE,TARGET :: dltxc(:)
   REAL(sp),ALLOCATABLE,TARGET :: dltyc(:)
   REAL(sp),ALLOCATABLE,TARGET :: dltxyc(:)
   REAL(sp),ALLOCATABLE,TARGET :: sitae(:) 
   REAL(sp),ALLOCATABLE,TARGET :: xijc(:) 
   REAL(sp),ALLOCATABLE,TARGET :: yijc(:)
   ! POSITION OF NODAL CONTROL VOLUME CORNERS 
   REAL(sp),ALLOCATABLE,TARGET :: xije(:,:) 
   REAL(sp),ALLOCATABLE,TARGET :: yije(:,:) 
   ! LENGTH OF NODAL CONTROL VOLUME EDGES
   REAL(sp),ALLOCATABLE,TARGET :: dltxe(:)
   REAL(sp),ALLOCATABLE,TARGET :: dltye(:)
   REAL(sp),ALLOCATABLE,TARGET :: dltxye(:)
   REAL(sp),ALLOCATABLE,TARGET :: sitac(:) 


   REAL(sp),ALLOCATABLE,TARGET :: epor(:)            !!ELEMENT FLUX POROSITY (=0. IF ISBCE = 2)

   ! LENGTH BETWEEN NODE AND CONTROL VOLUMEN EDGE CENTER
   REAL(sp),ALLOCATABLE,TARGET :: dltxncve(:,:)!! DeLTa X Node to Control Volume Edge 
   REAL(sp),ALLOCATABLE,TARGET :: dltyncve(:,:)!! DeLTa Y Node to Control Volume Edge 


   REAL(sp),ALLOCATABLE,TARGET :: dltytrie(:,:)!! DeLTa Y TRIangle Edge
   REAL(sp),ALLOCATABLE,TARGET :: dltxtrie(:,:)!! DeLTa X TRIangle Edge
   
   REAL(sp),ALLOCATABLE,TARGET :: dltxecec(:,:)!! DeLTa X Edge Center to Edge Center
   REAL(sp),ALLOCATABLE,TARGET :: dltyecec(:,:)!! DeLTa Y Edge Center to Edge Center

   REAL(sp),ALLOCATABLE,TARGET :: dltxnec(:,:)!! DeLTa X Node to Edge Center
   REAL(sp),ALLOCATABLE,TARGET :: dltynec(:,:)!! DeLTa Y Node to Edge Center


   ! LONG SHORE FLOW VARIABLES
   INTEGER, ALLOCATABLE,TARGET :: ibclsf_gl(:)     !!GLOBAL NODE NUMBER OF LSF BOUNDARY   
   REAL(sp),ALLOCATABLE,TARGET :: rbc_geo_gl(:)     !!GLOBAL GEOSTROPHIC FRICTION CORRECTION NODES
   REAL(sp),ALLOCATABLE,TARGET :: rbc_wdf_gl(:)     !!GLOBAL WIND DRIVEN FLOW CORRECTION NODES

   REAL(sp),ALLOCATABLE,TARGET :: wdf_ang(:)       !!ANGLE ALLONG THE OPEN BOUNDARY
   REAL(sp),ALLOCATABLE,TARGET :: wdf_dist(:)      !!DISTANCE ALLONG THE OPEN BOUNDARY 
   INTEGER, ALLOCATABLE,TARGET :: ibclsf(:)        !!LOCAL NODE NUMBER OF LSF BOUNDARY   
   INTEGER, ALLOCATABLE,TARGET :: ibclsf_output(:) !! LIST OF LOCAL LSF NODES GLOBAL NUMBER FOR OUTPUT   
   INTEGER, ALLOCATABLE,TARGET :: nbclsf(:)        !!LOCAL NODE NUMBER OF THE NEXT LSF BOUNDARY   
   REAL(sp),ALLOCATABLE,TARGET :: rbc_geo(:)        !!LOCAL GEOSTROPHIC FRICTION CORRECTION NODES
   REAL(sp),ALLOCATABLE,TARGET :: rbc_wdf(:)        !!LOCAL WIND DRIVEN FLOW CORRECTION NODES

!   INTEGER, ALLOCATABLE,TARGET :: N_ICELLQ(:,:)    !!FLUX ANGLE 

!----------------2-d arrays for the general vertical coordinate -------------------------------!

   REAL(sp), ALLOCATABLE,TARGET :: z(:,:)                    !!SIGMA COORDINATE VALUE 
   REAL(sp), ALLOCATABLE,TARGET :: zz(:,:)                   !!INTRA LEVEL SIGMA VALUE
   REAL(sp), ALLOCATABLE,TARGET :: dz(:,:)                   !!DELTA-SIGMA VALUE
   REAL(sp), ALLOCATABLE,TARGET :: dzz(:,:)                  !!DELTA OF INTRA LEVEL SIGMA 
   REAL(sp), ALLOCATABLE,TARGET :: z1(:,:)                   !!SIGMA COORDINATE VALUE 
   REAL(sp), ALLOCATABLE,TARGET :: zz1(:,:)                  !!INTRA LEVEL SIGMA VALUE
   REAL(sp), ALLOCATABLE,TARGET :: dz1(:,:)                  !!DELTA-SIGMA VALUE
   REAL(sp), ALLOCATABLE,TARGET :: dzz1(:,:)                 !!DELTA OF INTRA LEVEL SIGMA 
!   REAL(SP), ALLOCATABLE,TARGET :: DPTHSL(:)               !!Z-DEPTHS FOR SALINITY/TEMP ICs


!---------------2-d flow variable arrays at elements-------------------------------!

   REAL(sp), ALLOCATABLE,TARGET :: ua(:)            !!VERTICALLY AVERAGED X-VELOC
   REAL(sp), ALLOCATABLE,TARGET :: va(:)            !!VERTICALLY AVERAGED Y-VELOC
   REAL(sp), ALLOCATABLE,TARGET :: uaf(:)           !!UA FROM PREVIOUS RK STAGE 
   REAL(sp), ALLOCATABLE,TARGET :: vaf(:)           !!VA FROM PREVIOUS RK STAGE 
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: uark(:)          !!UA FROM PREVIOUS TIMESTEP 
   REAL(sp), ALLOCATABLE,TARGET :: vark(:)          !!VA FROM PREVIOUS TIMESTEP 
   REAL(sp), ALLOCATABLE,TARGET :: uard(:)          !!UA AVERAGED OVER EXTERNAL INT
   REAL(sp), ALLOCATABLE,TARGET :: vard(:)          !!VA AVERAGED OVER EXTERNAL INT
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: cor(:)           !!CORIOLIS PARAMETER
   REAL(sp), ALLOCATABLE,TARGET :: f_alfa(:)
   REAL(sp), ALLOCATABLE,TARGET :: h1(:)            !!BATHYMETRIC DEPTH   
   REAL(sp), ALLOCATABLE,TARGET :: d1(:)            !!CURRENT DEPTH
   REAL(sp), ALLOCATABLE,TARGET :: dt1(:)           !!DEPTH AT PREVIOUS TIME STEP
   REAL(sp), ALLOCATABLE,TARGET :: el1(:)           !!CURRENT SURFACE ELEVATION
   REAL(sp), ALLOCATABLE,TARGET :: et1(:)           !!SURFACE ELEVATION AT PREVIOUS TIME STEP
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: elrk1(:)         !!SURFACE ELEVATION AT BEGINNING OF RK INT 
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: elf1(:)          !!SURFACE ELEVATION STORAGE FOR RK INT
   REAL(sp), ALLOCATABLE,TARGET :: dtfa(:)          !!ADJUSTED DEPTH FOR MASS CONSERVATION


   REAL(sp), ALLOCATABLE,TARGET :: cc_sponge(:)     !!SPONGE DAMPING COEFFICIENT FOR MOMENTUM
   
!---------------2-d flow variable arrays at nodes----------------------------------!

   REAL(sp), ALLOCATABLE,TARGET :: h(:)     !!BATHYMETRIC DEPTH   
   REAL(sp), ALLOCATABLE,TARGET :: d(:)             !!CURRENT DEPTH   
   REAL(sp), ALLOCATABLE,TARGET :: dt(:)            !!DEPTH AT PREVIOUS TIME STEP
   REAL(sp), ALLOCATABLE,TARGET :: el(:)            !!CURRENT SURFACE ELEVATION
   REAL(sp), ALLOCATABLE,TARGET :: et(:)            !!SURFACE ELEVATION AT PREVIOUS TIME STEP
   REAL(sp), ALLOCATABLE,TARGET :: egf(:)           !!AVERAGE SURFACE ELEVATION OVER EXTERNAL INT
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: elrk(:)          !!SURFACE ELEVATION AT BEGINNING OF RK INT
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: elf(:)           !!SURFACE ELEVATION STORAGE FOR RK INT

   ! DEFINED HERE, BUT ONLY USED IF EQI&ATMO ARE DEFINED
   REAL(sp), ALLOCATABLE,TARGET :: elf_eqi(:)
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: elrk_eqi(:)
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: el_eqi(:)
   REAL(sp), ALLOCATABLE,TARGET :: egf_eqi(:)

   REAL(sp), ALLOCATABLE,TARGET :: elf_atmo(:)
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: elrk_atmo(:)
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: el_atmo(:)
   REAL(sp), ALLOCATABLE,TARGET :: egf_atmo(:)

   ! DEFINED HERE, BUT ONLY USED IF AIR PRESSURE ARE DEFINED
   REAL(sp), ALLOCATABLE,TARGET :: elf_air(:)
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: elrk_air(:)
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: el_air(:)
   REAL(sp), ALLOCATABLE,TARGET :: egf_air(:)
   
   REAL(sp), ALLOCATABLE,TARGET :: vort(:)


!---------------surface/bottom boundary conditions---------------------------------!

   REAL(sp), ALLOCATABLE,TARGET :: cbc(:)           !!BOTTOM FRICTION     
   REAL(sp), ALLOCATABLE,TARGET :: cc_z0b(:)        !!BOTTOM ROUGHNESS VARIABLE

   REAL(sp), ALLOCATABLE,TARGET :: swrad_watts(:)  !!SURFACE INCIDENT RADIATION
   REAL(sp), ALLOCATABLE,TARGET :: wtsurf_watts(:) !!NET HEAT FLUX AT SURFACE

   ! THESE TWO ARE DIVIDED BY THE SPECFIC HEAT AND AVERAGE DENSITY!
   REAL(sp), ALLOCATABLE,TARGET :: swrad(:)         !!SURFACE INCIDENT RADIATION
   REAL(sp), ALLOCATABLE,TARGET :: wtsurf(:)        !!NET HEAT FLUX AT SURFACE

   ! SUSPECTED UNITS FOR SURFACE STRESS: N/M * 1/RHO_BAR ... see bcond_gcn.F
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: wusurf2(:)       !!SURFACE FRICTION FOR EXT
   REAL(sp), ALLOCATABLE,TARGET :: wvsurf2(:)       !!SURFACE FRICTION FOR EXT
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: wubot(:)         !!BOTTOM FRICTION
   REAL(sp), ALLOCATABLE,TARGET :: wvbot(:)         !!BOTTOM FRICTION
   REAL(sp), ALLOCATABLE,TARGET :: taubm(:)         !!BOTTOM FRICTION MAGNITUDE(Caution is Tau' [no Rho])
   REAL(sp), ALLOCATABLE,TARGET :: wubot_n(:)       !!BOTTOM FRICTION ON NODES (Caution is Tau' [no Rho])
   REAL(sp), ALLOCATABLE,TARGET :: wvbot_n(:)       !!BOTTOM FRICTION ON NODES (Caution is Tau' [no Rho])
   REAL(sp), ALLOCATABLE,TARGET :: taubm_n(:)       !!BOTTOM FRICTION MAGNITUDE ON NODES (Caution is Tau' [no Rho])
   REAL(sp), ALLOCATABLE,TARGET :: wusurf(:)        !!SURFACE FRICTION FOR INT
   REAL(sp), ALLOCATABLE,TARGET :: wvsurf(:)        !!SURFACE FRICTION FOR INT
   REAL(sp), ALLOCATABLE,TARGET :: wusurf_save(:)   !!SURFACE FRICTION FOR INT
   REAL(sp), ALLOCATABLE,TARGET :: wvsurf_save(:)   !!SURFACE FRICTION FOR INT
   ! BFWDIS - UNITS: m3s-1 Cubic meters per second
   REAL(sp), ALLOCATABLE,TARGET :: bfwdis(:)        !!GROUNDWATER FLUX AT CURRENT TIME
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: bfwdis2(:)       !!GROUNDWATER FLUX FOR EXT
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: bfwtmp(:)        !!GROUNDWATER TEMP AT CURRENT TIME
   REAL(sp), ALLOCATABLE,TARGET :: bfwslt(:)        !!GROUNDWATER SALT AT CURRENT TIME


   ! ICING MODEL DATA
   REAL(sp), ALLOCATABLE,TARGET :: icing_wndx(:)
   REAL(sp), ALLOCATABLE,TARGET :: icing_wndy(:)
   REAL(sp), ALLOCATABLE,TARGET :: icing_satmp(:)
   REAL(sp), ALLOCATABLE,TARGET :: icing_0kts(:)
   REAL(sp), ALLOCATABLE,TARGET :: icing_10kts(:)

   ! NOT SURE WHETHER RIVER STUFF BELONGS HERE OR IN ALL_VARS???

   ! RIVER STUFF
   INTEGER,  ALLOCATABLE,TARGET :: inodeq(:)        !!LOCAL FRESH WATER INFLOW NODES
   INTEGER,  ALLOCATABLE,TARGET :: icellq(:)        !!LOCAL FRESH WATER INFLOW ELEMENT 
   INTEGER,  ALLOCATABLE,TARGET :: n_icellq(:,:)    !!TWO NODES BOUNDING THE INFLOW ELEMENT
   REAL(sp), ALLOCATABLE,TARGET :: vqdist(:,:)      !!DISCHARGE VERTICAL DISTRIBUTION
   INTEGER,  ALLOCATABLE,TARGET :: riv_gl2loc(:)

   REAL(sp), ALLOCATABLE,TARGET :: qdis(:)          !!RIVER FLUX AT CURRENT TIME
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: qdis2(:)         !!RIVER FLUX (EXT MODE, NOT USED) 
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: tdis(:)          !!RIVER WATER TEMP AT CURRENT TIME
   REAL(sp), ALLOCATABLE,TARGET :: sdis(:)          !!RIVER WATER SLNT AT CURRENT TIME
   REAL(sp), ALLOCATABLE,TARGET :: qarea(:)         !!AREA OF RIVER DISCHARGE 
   REAL(sp), ALLOCATABLE,TARGET :: rdisq(:,:)       !!AREA OF FLUX                    
   REAL(sp), ALLOCATABLE,TARGET :: angleq(:)        !!RIVER DISCHARGE ANGLE           
   REAL(sp), ALLOCATABLE,TARGET :: vlctyq(:)        !!RIVER DISCHARGE VELOCITY
   
   ! SURFACE MET STUFF
   REAL(sp), ALLOCATABLE,TARGET :: uuwind(:)        !!SURFACE X-WIND 
   REAL(sp), ALLOCATABLE,TARGET :: vvwind(:)        !!SURFACE Y-WIND
   ! PRECIP/EVAP are in units of meters/second
!# if !defined (SEMI_IMPLICIT)
   REAL(sp), ALLOCATABLE,TARGET :: qprec2(:)        !!SURFACE PRECIPITATION FOR EXT
   REAL(sp), ALLOCATABLE,TARGET :: qevap2(:)        !!SURFACE EVAPORATION FOR EXT
!# endif
   REAL(sp), ALLOCATABLE,TARGET :: qprec(:)        !!SURFACE PRECIPITATION FOR INT
   REAL(sp), ALLOCATABLE,TARGET :: qevap(:)        !!SURFACE EVAPORATION FOR INT

   REAL(sp), ALLOCATABLE,TARGET :: whs(:)          !!SURFACE WAVE HEIGHT
   REAL(sp), ALLOCATABLE,TARGET :: wdir(:)         !!SURFACE WAVE DIRECTION
   REAL(sp), ALLOCATABLE,TARGET :: wper(:)         !!SURFACE WAVE PERIOD
   REAL(sp), ALLOCATABLE,TARGET :: wlength(:)      !!SURFACE WAVE LENGTH
   REAL(sp), ALLOCATABLE,TARGET :: wper_bot(:)     !!BOTTOM WAVE PERIOD
   REAL(sp), ALLOCATABLE,TARGET :: wub_bot(:)      !!BOTTOM ORBITAL VELOCITY


!----------------boundary conditions: meteo conditions-----------------!

!-----------------------2-d flow fluxes--------------------------------------------!

   REAL(sp), ALLOCATABLE,TARGET :: pstx(:)           !!EXT MODE BAROTROPIC TERMS
   REAL(sp), ALLOCATABLE,TARGET :: psty(:)           !!EXT MODE BAROTROPIC TERMS
   REAL(sp), ALLOCATABLE,TARGET :: advua(:) 
   REAL(sp), ALLOCATABLE,TARGET :: advva(:) 
   REAL(sp), ALLOCATABLE,TARGET :: adx2d(:) 
   REAL(sp), ALLOCATABLE,TARGET :: ady2d(:) 
   REAL(sp), ALLOCATABLE,TARGET :: drx2d(:) 
   REAL(sp), ALLOCATABLE,TARGET :: dry2d(:) 
   REAL(sp), ALLOCATABLE,TARGET :: tps(:)            !!WORKING ARRAY             
   REAL(sp), ALLOCATABLE,TARGET :: advx(:,:)      
   REAL(sp), ALLOCATABLE,TARGET :: advy(:,:)     

!---------------- internal mode   arrays-(element based)----------------------------!

   REAL(sp), ALLOCATABLE,TARGET :: u(:,:)         !X-VELOCITY
   REAL(sp), ALLOCATABLE,TARGET :: v(:,:)         !Y-VELOCITY

   REAL(sp), ALLOCATABLE,TARGET :: ubeta(:,:)     !X-VELOCITY temp time step
   REAL(sp), ALLOCATABLE,TARGET :: vbeta(:,:)     !Y-VELOCITY temp time step

   REAL(sp), ALLOCATABLE :: ubeta2d(:)
   REAL(sp), ALLOCATABLE :: vbeta2d(:)

   REAL(sp), ALLOCATABLE, TARGET :: partition(:) !gwc

   REAL(sp), ALLOCATABLE,TARGET :: w(:,:)         !VERTICAL VELOCITY IN SIGMA SYSTEM
   REAL(sp), ALLOCATABLE,TARGET :: ww(:,:)        !Z-VELOCITY
   REAL(sp), ALLOCATABLE,TARGET :: uf(:,:)        !X-VELOCITY FROM PREVIOUS TIMESTEP
   REAL(sp), ALLOCATABLE,TARGET :: vf(:,:)        !Y-VELOCITY FROM PREVIOUS TIMESTEP
   REAL(sp), ALLOCATABLE,TARGET :: wt(:,:)        !Z-VELOCITY FROM PREVIOUS TIMESTEP
   REAL(sp), ALLOCATABLE,TARGET :: rho(:,:)       !DENSITY AT ELEMENTS
   REAL(sp), ALLOCATABLE,TARGET :: rmean(:,:)     !INITIAL DENSITY AT ELEMENTS
   REAL(sp), ALLOCATABLE,TARGET :: t(:,:)         !TEMPERATURE AT ELEMENTS
   REAL(sp), ALLOCATABLE,TARGET :: tmean(:,:)     !INITIAL TEMPERATURE AT ELEMENTS
   REAL(sp), ALLOCATABLE,TARGET :: s(:,:)         !SALINITY AT ELEMENTS
   REAL(sp), ALLOCATABLE,TARGET :: smean(:,:)     !INITIAL SALINITY AT ELEMENTS
   REAL(sp), ALLOCATABLE,TARGET :: q2(:,:)        !2 X TURBULENT KINETIC ENERGY AT NODES
   REAL(sp), ALLOCATABLE,TARGET :: l(:,:)         !TURBULENT LENGTH MACROSCALE 
   REAL(sp), ALLOCATABLE,TARGET :: q2l(:,:)       !2 X TURBULENT KE X LENGTH AT NODES
  REAL(sp), ALLOCATABLE,TARGET :: km(:,:)        !TURBULENT EDDY VISCOSITY FOR MOMENTUM 
  REAL(sp), ALLOCATABLE,TARGET :: kh(:,:)        !TURBULENT DIFFUSIVITY FOR SALINITY/TEMP 
  REAL(sp), ALLOCATABLE,TARGET :: kq(:,:)        !TURBULENT DIFFUSIVITY FOR Q2/Q2L 
  REAL(sp), ALLOCATABLE,TARGET :: aam(:,:)       !STORAGE FOR OUTPUT OF HORIZONTAL VISCOSITY
  REAL(sp), ALLOCATABLE,TARGET :: q2f(:,:)       !WORKING ARRAY FOR UPDATING Q2
  REAL(sp), ALLOCATABLE,TARGET :: q2lf(:,:)      !WORKING ARRAY FOR UPDATING Q2F 
  REAL(sp), ALLOCATABLE,TARGET :: km1(:,:)       !TURBULENT EDDY VISCOSITY FOR MOMENTUM 

  ! VARIABLE HORIZONTAL VISCOSITY COEFFICENTS  
  REAL(sp), ALLOCATABLE :: cc_hvc(:)
  REAL(sp), ALLOCATABLE :: nn_hvc(:)

  !-----------------------3d variable arrays-(node based)-----------------------------!

  REAL(sp), ALLOCATABLE,TARGET :: t1(:,:)         !!TEMPERATURE AT NODES               
  REAL(sp), ALLOCATABLE,TARGET :: s1(:,:)         !!SALINITY AT NODES               
  REAL(sp), ALLOCATABLE,TARGET :: rho1(:,:)       !!DENSITY AT NODES               
  REAL(sp), ALLOCATABLE,TARGET :: tf1(:,:)        !!TEMPERATURE FROM PREVIOUS TIME
  REAL(sp), ALLOCATABLE,TARGET :: sf1(:,:)        !!SALINITY FROM PREVIOUS TIME 
!J. Ge for tracer advection
  REAL(sp), ALLOCATABLE,TARGET :: t0(:,:)         !!TEMPERATURE AT NODES AT PREVIOUS TIME STEP
  REAL(sp), ALLOCATABLE,TARGET :: t2(:,:)         !!TEMPERATURE AT NODES AT PREVIOUS TWO STEP
  REAL(sp), ALLOCATABLE,TARGET :: s0(:,:)         !!SALINITY AT NODES AT PREVIOUS TIME STEP
  REAL(sp), ALLOCATABLE,TARGET :: s2(:,:)         !!SALINITY AT NODES AT PREVIOUS TWO STEP
!J. Ge for tracer advection
  REAL(sp), ALLOCATABLE,TARGET :: tmean1(:,:)     !!MEAN INITIAL TEMP
  REAL(sp), ALLOCATABLE,TARGET :: smean1(:,:)     !!MEAN INITIAL SALINITY 
  REAL(sp), ALLOCATABLE,TARGET :: rmean1(:,:)     !!MEAN INITIAL DENSITY 
  REAL(sp), ALLOCATABLE,TARGET :: wts(:,:)        !!VERTICAL VELOCITY IN SIGMA SYSTEM
  REAL(sp), ALLOCATABLE,TARGET :: wtts(:,:)       !!WTS FROM PREVIOUS TIMESTEP        

  !---------------------------internal mode fluxes------------------------------------!

  REAL(sp), ALLOCATABLE,TARGET :: drhox(:,:)      !!BAROCLINIC PG IN X DIRECTION
  REAL(sp), ALLOCATABLE,TARGET :: drhoy(:,:)      !!BAROCLINIC PG IN Y DIRECTION

  !------------shape coefficient arrays and control volume metrics--------------------!

  REAL(sp), ALLOCATABLE,TARGET :: a1u(:,:)      
  REAL(sp), ALLOCATABLE,TARGET :: a2u(:,:)     
  REAL(sp), ALLOCATABLE,TARGET :: awx(:,:)   
  REAL(sp), ALLOCATABLE,TARGET :: awy(:,:)  
  REAL(sp), ALLOCATABLE,TARGET :: aw0(:,:) 
  REAL(sp), ALLOCATABLE,TARGET :: alpha(:)


  !-----salinity and temperature bottom diffusion condition/bottom depth gradients----!

  REAL(sp), ALLOCATABLE,TARGET :: phpn(:)   
  REAL(sp), ALLOCATABLE,TARGET :: pfpxb(:)   
  REAL(sp), ALLOCATABLE,TARGET :: pfpyb(:)   
  REAL(sp), ALLOCATABLE,TARGET :: sita_gd(:)  
  REAL(sp), ALLOCATABLE,TARGET :: ah_bottom(:)  

  !-----arrays used for averaging of flow quantities for output-----------------------!

  REAL(sp), ALLOCATABLE,TARGET :: u_ave(:,:)   !U AVERAGED OVER INT_AVGE ITERATIONS
  REAL(sp), ALLOCATABLE,TARGET :: v_ave(:,:)   !V AVERAGED OVER INT_AVGE ITERATIONS
  REAL(sp), ALLOCATABLE,TARGET :: w_ave(:,:)   !WW AVERAGED OVER INT_AVGE ITERATIONS
  REAL(sp), ALLOCATABLE,TARGET :: km_ave(:,:)  !KM AVERAGED OVER INT_AVGE ITERATIONS
  REAL(sp), ALLOCATABLE,TARGET :: kh_ave(:,:)  !KH AVERAGED OVER INT_AVGE ITERATIONS
  REAL(sp), ALLOCATABLE,TARGET :: t_ave(:,:)   !T1 AVERAGED OVER INT_AVGE ITERATIONS
  REAL(sp), ALLOCATABLE,TARGET :: s_ave(:,:)   !S1 AVERAGED OVER INT_AVGE ITERATIONS
  REAL(sp), ALLOCATABLE,TARGET :: r_ave(:,:)   !RHO1 AVERAGED OVER INT_AVGE ITERATIONS
  REAL(sp), ALLOCATABLE,TARGET :: el_ave(:)    !EL AVERAGED OVER INT_AVGE ITERATIONS

  REAL(sp), ALLOCATABLE,TARGET :: viscofh(:,:)
  REAL(sp), ALLOCATABLE,TARGET :: viscofm(:,:)

  REAL(sp), ALLOCATABLE,TARGET :: hyw(:,:)


CONTAINS


!==============================================================================|
   SUBROUTINE n2e3d(NVAR,EVAR)
!==============================================================================|
   IMPLICIT NONE 
   REAL(SP), DIMENSION(0:MT,1:KB), INTENT(IN)  :: NVAR  
   REAL(SP), DIMENSION(0:NT,1:KB), INTENT(INOUT) :: EVAR  
   INTEGER I,K
!------------------------------------------------------------------------------|
   DO k=1,kb
     DO i = 1, nt
       evar(i,k) = one_third*(nvar(nv(i,1),k)+nvar(nv(i,2),k)+nvar(nv(i,3),k))
     END DO
   END DO
   RETURN
   END SUBROUTINE n2e3d 
!==============================================================================|

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!==============================================================================|
   SUBROUTINE n2e2d(NVAR,EVAR)
!==============================================================================|
   IMPLICIT NONE
   REAL(SP), DIMENSION(0:MT), INTENT(IN)  :: NVAR 
   REAL(SP), DIMENSION(0:NT), INTENT(INOUT) :: EVAR  
   INTEGER I,K
!------------------------------------------------------------------------------|
   DO i = 1, nt
     evar(i) = one_third*(nvar(nv(i,1))+nvar(nv(i,2))+nvar(nv(i,3)))
   END DO
   RETURN
   END SUBROUTINE n2e2d
!==============================================================================|

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!==============================================================================|
   SUBROUTINE e2n2d(EVAR,NVAR)
!==============================================================================|
   IMPLICIT NONE
   REAL(SP), DIMENSION(0:NT), INTENT(IN ) :: EVAR  
   REAL(SP), DIMENSION(0:MT), INTENT(INOUT) :: NVAR 

   INTEGER I,K
!------------------------------------------------------------------------------|
   DO i=1,m
     nvar(i) = sum(evar(nbve(i,1:ntve(i))))/float(ntve(i))
   END DO
   RETURN
   END SUBROUTINE e2n2d
!==============================================================================|

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!==============================================================================|
   SUBROUTINE e2n3d(EVAR,NVAR)
!==============================================================================|
   IMPLICIT NONE
   REAL(SP), DIMENSION(0:NT,1:KB), INTENT(IN ) :: EVAR  
   REAL(SP), DIMENSION(0:MT,1:KB), INTENT(INOUT) :: NVAR 
   INTEGER I,K
!------------------------------------------------------------------------------|
   DO k=1,kb
     DO i=1,m
       nvar(i,k) = sum(evar(nbve(i,1:ntve(i)),k))/float(ntve(i))
     END DO
   END DO
   RETURN
   END SUBROUTINE e2n3d
!==============================================================================|
   
!==============================================================================|
!    Allocate and Initialize Most Arrays                                       !
!==============================================================================|
   
   SUBROUTINE alloc_vars(dbg_set)
     
!==============================================================================!
   IMPLICIT NONE
   logical,intent(in) :: dbg_set
   INTEGER NCT,NDB
!==============================================================================!
   ndb = 1       !!GWC BASE THIS ON KIND

   nct = nt*3 ! A DIMENSION USED FOR ALLOCATION!
   
!==============================================================================!
!  ALLOCATE:                                                                   !
!==============================================================================!
   
!--------------------------Grid Metrics---------------------------------------------!
   
   ALLOCATE(lon(0:mt))           ;lon  = zero   !!LONGITUDE AT THE NODE
   ALLOCATE(lat(0:mt))           ;lat  = zero   !!LATITUDE AT THE NODE
   ALLOCATE(lonc(0:nt))          ;lonc = zero   !!LONGITUDE AT THE NODE
   ALLOCATE(latc(0:nt))          ;latc = zero   !!LATITUDE AT THE NODE
   ALLOCATE(xm(0:mt))            ;xm   = zero   !!X-COORD AT NODE IN METERS
   ALLOCATE(ym(0:mt))            ;ym   = zero   !!Y-COORD AT NODE IN METERS
   ALLOCATE(xmc(0:nt))           ;xmc  = zero   !!X-COORD AT FACE CENTER IN METERS 
   ALLOCATE(ymc(0:nt))           ;ymc  = zero   !!Y-COORD AT FACE CENTER IN METERS
   ALLOCATE(xc(0:nt))            ;xc   = zero   !!X-COORD AT FACE CENTER 
   ALLOCATE(yc(0:nt))            ;yc   = zero   !!Y-COORD AT FACE CENTER
   ALLOCATE(vx(0:mt))            ;vx   = zero   !!X-COORD AT GRID POINT
   ALLOCATE(vy(0:mt))            ;vy   = zero   !!Y-COORD AT GRID POINT
   ALLOCATE(art(0:nt))           ;art  = zero   !!AREA OF ELEMENT
   ALLOCATE(art1(0:mt))          ;art1 = zero   !!AREA OF NODE-BASE CONTROl VOLUME
   ALLOCATE(art2(0:mt))          ;art2 = zero   !!AREA OF ELEMENTS AROUND NODE
   ALLOCATE(grav_n(0:mt))        ;grav_n = zero !! VARIABLE GRAVITY
   ALLOCATE(grav_e(0:nt))        ;grav_e = zero !! VARIABLE GRAVITY

   memcnt = mt*9*ndb + nt*8*ndb + memcnt

!----------------Node, Boundary Condition, and Control Volume-----------------------!

! ALLOCATED IN setup_domain.F But count the memory here
!   ALLOCATE(NV(0:NT,4))          ;NV       = 0  !!NODE NUMBERING FOR ELEMENTS
!   ALLOCATE(NVGL(0:NT,3))        ;NVGL     = 0  !!GLOBAL NODE NUMBERING FOR LOCAL ELEMENTS
   ALLOCATE(nbe(0:nt,3))         ;nbe      = 0 !!INDICES OF ELMNT NEIGHBORS
   ALLOCATE(ntve(0:mt))          ;ntve     = 0 
   ALLOCATE(ntsn(0:mt))          ;ntsn     = 0 
   ALLOCATE(isonb(0:mt))         ;isonb    = 0  !!NODE MARKER = 0,1,2
   ALLOCATE(isonb_w(0:mt))       ;isonb_w  = 0  !!NODE MARKER = 0,1,2
   ALLOCATE(isbce(0:nt))         ;isbce    = 0 
   ALLOCATE(niec(nct,2))         ;niec     = 0
   ALLOCATE(ntrg(nct))           ;ntrg     = 0

   ! POSITION OF NODAL CONTROL VOLUME CORNERS 
   ALLOCATE(xije(nct,2))         ;xije     = zero
   ALLOCATE(yije(nct,2))         ;yije     = zero 
   ! LENGTH OF NODAL CONTROL VOLUME EDGES
   ALLOCATE(dltxe(nct))          ;dltxe    = zero
   ALLOCATE(dltye(nct))          ;dltye    = zero
   ALLOCATE(dltxye(nct))         ;dltxye   = zero !! TOTAL LENGTH
   ALLOCATE(sitae(nct))          ;sitae    = zero !! ANGLE

   ! LENGTH BETWEEN NODE AND CONTROL VOLUMEN EDGE CENTER
   ALLOCATE(dltxncve(nct,2))     ; dltxncve   = zero  !! DeLTa X Node to Control Volume Edge 
   ALLOCATE(dltyncve(nct,2))     ; dltyncve   = zero  !! DeLTa Y Node to Control Volume Edge 


   ! THE FOLLOWING ARRAYS COULD BE REPLACED WITH (N,3) ARRAYS
   ! BUT THE INDEXING WOULD BE COMPLEX AND THERE ARE SIGN ISSUES!

   ! TRIANGLE EDGE LENGTH FOR EDGES SURROUNDING EACH NODE
   ! NTSN Can not be greater than 13!
   ALLOCATE(dltxtrie(m,12))         ;dltxtrie   = zero !! DeLTa X TRIangle Edge
   ALLOCATE(dltytrie(m,12))         ;dltytrie   = zero !! DeLTa Y TRIangle Edge

   ALLOCATE(dltxnec(m,12))         ;dltxnec   = zero !! DeLTa X Node to Edge Center
   ALLOCATE(dltynec(m,12))         ;dltynec   = zero !! DeLTa Y Node to Edge Center

   ! DISTANCE BETWEEN TRIANGLE EDGE CENTERS FOR EACH TRIANGLE AROUND A NODE
   ! NTVE Can not be greater than 13!  
   ALLOCATE(dltxecec(m,12))         ;dltxecec   = zero  !! DeLTa X Edge Center to Edge Center
   ALLOCATE(dltyecec(m,12))         ;dltyecec   = zero  !! DeLTa Y Edge Center to Edge Center


   
   memcnt = nt*4 + mt*3 +m*6*12*ndb  + nct*3 + nct*12*ndb  + memcnt

!----------------2-d arrays for the general vertical coordinate -------------------------------!

   ALLOCATE(z(0:mt,kb))               ; z      = zero    !!SIGMA COORDINATE VALUE 
   ALLOCATE(zz(0:mt,kb))              ; zz     = zero    !!INTRA LEVEL SIGMA VALUE
   ALLOCATE(dz(0:mt,kb))              ; dz     = zero    !!DELTA-SIGMA VALUE
   ALLOCATE(dzz(0:mt,kb))             ; dzz    = zero    !!DELTA OF INTRA LEVEL SIGMA 
   ALLOCATE(z1(0:nt,kb))              ; z1     = zero    !!SIGMA COORDINATE VALUE 
   ALLOCATE(zz1(0:nt,kb))             ; zz1    = zero    !!INTRA LEVEL SIGMA VALUE
   ALLOCATE(dz1(0:nt,kb))             ; dz1    = zero    !!DELTA-SIGMA VALUE
   ALLOCATE(dzz1(0:nt,kb))            ; dzz1   = zero    !!DELTA OF INTRA LEVEL SIGMA 
   memcnt = mt*kb*4*ndb + nt*kb*4*ndb +memcnt

!---------------2-d flow variable arrays at elements-------------------------------!

   ALLOCATE(ua(0:nt))            ;ua        = zero  !!VERTICALLY AVERAGED X-VELOC
   ALLOCATE(va(0:nt))            ;va        = zero  !!VERTICALLY AVERAGED Y-VELOC
   ALLOCATE(uaf(0:nt))           ;uaf       = zero  !!UA FROM PREVIOUS RK STAGE 
   ALLOCATE(vaf(0:nt))           ;vaf       = zero  !!VA FROM PREVIOUS RK STAGE 
   ALLOCATE(uark(0:nt))          ;uark      = zero  !!UA FROM PREVIOUS TIMESTEP 
   ALLOCATE(vark(0:nt))          ;vark      = zero  !!VA FROM PREVIOUS TIMESTEP 
   ALLOCATE(uard(0:nt))          ;uard      = zero  !!UA AVERAGED OVER EXTERNAL INT
   ALLOCATE(vard(0:nt))          ;vard      = zero  !!VA AVERAGED OVER EXTERNAL INT
   ALLOCATE(cor(0:nt))           ;cor       = zero  !!CORIOLIS PARAMETER
   ALLOCATE(f_alfa(0:nt))        ;f_alfa    = 1.0_sp  !!EQUATORIAL BETA PLANE PARAMETER
   ALLOCATE(h1(0:nt))            ;h1        = zero  !!BATHYMETRIC DEPTH   
   ALLOCATE(d1(0:nt))            ;d1        = zero  !!DEPTH
   ALLOCATE(dt1(0:nt))           ;dt1       = zero  !!DEPTH  
   ALLOCATE(el1(0:nt))           ;el1       = zero  !!SURFACE ELEVATION
   ALLOCATE(elf1(0:nt))          ;elf1      = zero  !!SURFACE ELEVATION
   ALLOCATE(dtfa(0:mt))          ;dtfa      = zero  !!ADJUSTED DEPTH FOR MASS CONSERVATION
   ALLOCATE(et1(0:nt))           ;et1       = zero  !!SURFACE ELEVATION
   ALLOCATE(elrk1(0:nt))         ;elrk1     = zero  !!SURFACE ELEVATION
   ALLOCATE(cc_sponge(0:nt))     ;cc_sponge = zero  !!SPONGE DAMPING COEFFICIENT FOR MOMENTUM
                 memcnt = nt*17*ndb + mt*ndb + memcnt

!---------------2-d flow variable arrays at nodes----------------------------------!

   ALLOCATE(h(0:mt))             ;h    = zero       !!BATHYMETRIC DEPTH   
   ALLOCATE(d(0:mt))             ;d    = zero       !!DEPTH   
   ALLOCATE(dt(0:mt))            ;dt   = zero       !!DEPTH   
   ALLOCATE(el(0:mt))            ;el   = zero       !!SURFACE ELEVATION
   ALLOCATE(elf(0:mt))           ;elf  = zero       !!SURFACE ELEVATION
   ALLOCATE(et(0:mt))            ;et   = zero       !!SURFACE ELEVATION
   ALLOCATE(egf(0:mt))           ;egf  = zero       !!SURFACE ELEVATION
   ALLOCATE(elrk(0:mt))          ;elrk = zero       !!SURFACE ELEVATION
               memcnt = mt*8*ndb + memcnt


   

   ALLOCATE(vort(0:mt))          ; vort     = zero            
               memcnt = mt*ndb + memcnt

!---------------surface/bottom/edge boundary conditions-----------------------------!

   ALLOCATE(cbc(0:nt))           ;cbc    = zero     !!BOTTOM FRICTION     
   ALLOCATE(cc_z0b(0:nt))        ;cc_z0b = zero     !!BOTTOM ROUGHNESS VARIABLE     
!# if !defined (SEMI_IMPLICIT)
   ALLOCATE(wusurf2(0:nt))       ;wusurf2= zero     !!SURFACE FRICTION FOR EXT
   ALLOCATE(wvsurf2(0:nt))       ;wvsurf2= zero     !!SURFACE FRICTION FOR EXT
!# endif
   ALLOCATE(wubot(0:nt))         ;wubot  = zero     !!BOTTOM FRICTION
   ALLOCATE(wvbot(0:nt))         ;wvbot  = zero     !!BOTTOM FRICTION
   ALLOCATE(taubm(0:nt))         ;taubm  = zero     !!BOTTOM FRICTION
   ALLOCATE(wubot_n(0:mt))       ;wubot_n  = zero   !!U-Component bottom shear stress on nodes  
   ALLOCATE(wvbot_n(0:mt))       ;wvbot_n  = zero   !!V-Component bottom shear stress on nodes
   ALLOCATE(taubm_n(0:mt))       ;taubm_n  = zero   !!Magnitude bottom shear stress on nodes
   ALLOCATE(wusurf(0:nt))        ;wusurf = zero     !!SURFACE FRICTION FOR INT
   ALLOCATE(wvsurf(0:nt))        ;wvsurf = zero     !!SURFACE FRICTION FOR INT
   ALLOCATE(wusurf_save(0:nt))   ;wusurf_save = zero!!SURFACE FRICTION FOR INT
   ALLOCATE(wvsurf_save(0:nt))   ;wvsurf_save = zero!!SURFACE FRICTION FOR INT
   ALLOCATE(uuwind(0:nt))        ;uuwind = zero     !!SURFACE X-WIND 
   ALLOCATE(vvwind(0:nt))        ;vvwind = zero     !!SURFACE Y-WIND 
   ALLOCATE(swrad(0:mt))         ;swrad  = zero     !!SURFACE INCIDENT RADIATION
   ALLOCATE(wtsurf(0:mt))        ;wtsurf = zero 
   ALLOCATE(swrad_watts(0:mt))  ;swrad_watts  = zero !!SURFACE INCIDENT RADIATION
   ALLOCATE(wtsurf_watts(0:mt)) ;wtsurf_watts = zero 

   ALLOCATE(qprec2(0:mt))        ;qprec2 = zero     !!SURFACE PRECIPITATION FOR EXT
   ALLOCATE(qevap2(0:mt))        ;qevap2 = zero     !!SURFACE EVAPORATION FOR EXT
   ALLOCATE(qprec(0:mt))         ;qprec = zero     !!SURFACE PRECIPITATION FOR INT
   ALLOCATE(qevap(0:mt))         ;qevap = zero     !!SURFACE EVAPORATION FOR INT



   memcnt = nt*10*ndb + mt*8*ndb + memcnt

   IF (icing_model) THEN
      ALLOCATE(icing_wndx(0:mt)) ;icing_wndx = zero
      ALLOCATE(icing_wndy(0:mt)) ;icing_wndy = zero
      ALLOCATE(icing_satmp(0:mt));icing_satmp  = zero
      ALLOCATE(icing_0kts(0:mt)) ;icing_0kts = zero
      ALLOCATE(icing_10kts(0:mt));icing_10kts= zero

      memcnt = memcnt + mt*4*ndb
   END IF
      

!--------------------------------------------------------------
!--------------------------------------------------------------

   ALLOCATE(bfwdis(0:mt))        ;bfwdis = zero     !!GROUNDWATER FLUX FOR INT
   ALLOCATE(bfwdis2(0:mt))       ;bfwdis2= zero     !!GROUNDWATER FLUX FOR EXT
   ALLOCATE(bfwslt(0:mt))        ;bfwslt = zero     !!GROUNDWATER SALT AT CURRENT TIME
   ALLOCATE(bfwtmp(0:mt))        ;bfwtmp = zero     !!GROUNDWATER TEMP AT CURRENT TIME

   memcnt =  mt*4*ndb + memcnt

!-----------------------2-d flow fluxes---------------------------------------------!

   ALLOCATE(pstx(0:nt))          ;pstx  = zero       !!EXT MODE BAROTROPIC TERMS
   ALLOCATE(psty(0:nt))          ;psty  = zero       !!EXT MODE BAROTROPIC TERMS
   ALLOCATE(advua(0:nt))         ;advua = zero 
   ALLOCATE(advva(0:nt))         ;advva = zero
   ALLOCATE(adx2d(0:nt))         ;adx2d = zero
   ALLOCATE(ady2d(0:nt))         ;ady2d = zero
   ALLOCATE(drx2d(0:nt))         ;drx2d = zero
   ALLOCATE(dry2d(0:nt))         ;dry2d = zero
   ALLOCATE(advx(0:nt,kb))       ;advx  = zero 
   ALLOCATE(advy(0:nt,kb))       ;advy  = zero 
   ALLOCATE(tps(0:nt))           ;tps   = zero      !!WORKING ARRAY             
   memcnt = nt*9*ndb + nt*kb*2*ndb + memcnt


!---------------- internal mode   arrays-(element based)----------------------------!

   ALLOCATE(u(0:nt,kb))          ;u     = zero   !!X-VELOCITY
   ALLOCATE(v(0:nt,kb))          ;v     = zero   !!Y-VELOCITY
   ALLOCATE(ubeta(0:nt,kb))      ;ubeta = zero   !!X-VELOCITY temp time step
   ALLOCATE(vbeta(0:nt,kb))      ;vbeta = zero   !!X-VELOCITY temp time step
   ALLOCATE(ubeta2d(0:nt))       ;ubeta2d = zero
   ALLOCATE(vbeta2d(0:nt))       ;vbeta2d = zero

   ALLOCATE(w(0:nt,kb))          ;w     = zero   !!VERTICAL VELOCITY IN SIGMA SYSTEM
   ALLOCATE(ww(0:nt,kb))         ;ww    = zero   !!Z-VELOCITY
   ALLOCATE(uf(0:nt,kb))         ;uf    = zero   !!X-VELOCITY FROM PREVIOUS TIMESTEP
   ALLOCATE(vf(0:nt,kb))         ;vf    = zero   !!Y-VELOCITY FROM PREVIOUS TIMESTEP
   ALLOCATE(wt(0:nt,kb))         ;wt    = zero   !!Z-VELOCITY FROM PREVIOUS TIMESTEP
   ALLOCATE(rho(0:nt,kb))        ;rho   = zero   !!DENSITY AT ELEMENTS
   ALLOCATE(rmean(0:nt,kb))      ;rmean = zero   !!MEAN INITIAL DENSITY AT ELEMENTS
   ALLOCATE(t(0:nt,kb))          ;t     = zero   !!TEMPERATURE AT ELEMENTS
   ALLOCATE(tmean(0:nt,kb))      ;tmean = zero   !!MEAN INITIAL TEMPERATURE AT ELEMENTS
   ALLOCATE(s(0:nt,kb))          ;s     = zero   !!SALINITY AT ELEMENTS
   ALLOCATE(smean(0:nt,kb))      ;smean = zero   !!MEAN INITIAL SALINITY AT ELEMENTS
               memcnt = nt*kb*13*ndb + memcnt

!-----------------------3d variable arrays-(node based)-----------------------------!

   ALLOCATE(t1(0:mt,kb))         ;t1     = zero  !!TEMPERATURE AT NODES               
   ALLOCATE(s1(0:mt,kb))         ;s1     = zero  !!SALINITY AT NODES               
!J. Ge for tracer advection
   ALLOCATE(t0(0:mt,kb))         ;t0    = zero   !!TEMPERATURE FROM PREVIOUS TIME STEP 
   ALLOCATE(t2(0:mt,kb))         ;t2    = zero   !!TEMPERATURE FROM PREVIOUS TWO STEP
   ALLOCATE(s0(0:mt,kb))         ;s0    = zero   !!SALINITY FROM PREVIOUS TIME STEP 
   ALLOCATE(s2(0:mt,kb))         ;s2    = zero   !!SALINITY FROM PREVIOUS TWO STEP 
!J. Ge for tracer advection
   ALLOCATE(rho1(0:mt,kb))       ;rho1   = zero  !!DENSITY AT NODES               
   ALLOCATE(tf1(0:mt,kb))        ;tf1    = zero  !!TEMPERATURE FROM PREVIOUS TIME
   ALLOCATE(sf1(0:mt,kb))        ;sf1    = zero  !!SALINITY FROM PREVIOUS TIME 
   ALLOCATE(tmean1(0:mt,kb))     ;tmean1 = zero  !!MEAN INITIAL TEMP
   ALLOCATE(smean1(0:mt,kb))     ;smean1 = zero  !!MEAN INITIAL SALINITY 
   ALLOCATE(rmean1(0:mt,kb))     ;rmean1 = zero  !!MEAN INITIAL DENSITY 
   ALLOCATE(wts(0:mt,kb))        ;wts    = zero  !!VERTICAL VELOCITY IN SIGMA SYSTEM
   ALLOCATE(wtts(0:mt,kb))       ;wtts   = zero  !!WTS FROM PREVIOUS TIMESTEP        
   ALLOCATE(q2(0:mt,kb))         ;q2    = zero   !!TURBULENT KINETIC ENERGY AT NODES
   ALLOCATE(q2l(0:mt,kb))        ;q2l   = zero   !!TURBULENT KE*LENGTH AT NODES
   ALLOCATE(l(0:mt,kb))          ;l     = zero   !!TURBULENT LENGTH SCALE AT ELEMENTS
   ALLOCATE(km(0:mt,kb))         ;km    = zero   !!TURBULENT QUANTITY
   ALLOCATE(kh(0:mt,kb))         ;kh    = zero   !!TURBULENT QUANTITY
   ALLOCATE(kq(0:mt,kb))         ;kq    = zero   !!TURBULENT QUANTITY
   ALLOCATE(aam(0:mt,kb))        ;aam   = zero   !!??
   ALLOCATE(km1(0:nt,kb))        ;km1   = zero   !!TURBULENT QUANTITY AT ELEMENTS


   ALLOCATE(cc_hvc(0:nt))        ;cc_hvc = zero  !!VISCOSITY COEFFICIENT AT ELEMENTS
   ALLOCATE(nn_hvc(0:mt))        ;nn_hvc = zero  !!VISCOSITY COEFFICIENT AT NODES

   memcnt = mt*kb*18*ndb + nt*kb*ndb + memcnt + (mt+nt)*ndb

  
!---------------------------internal mode fluxes------------------------------------!

   ALLOCATE(drhox(0:nt,kb))      ;drhox  = zero 
   ALLOCATE(drhoy(0:nt,kb))      ;drhoy  = zero 
   ALLOCATE(q2f(0:mt,kb))        ;q2f    = zero 
   ALLOCATE(q2lf(0:mt,kb))       ;q2lf   = zero
   memcnt = nt*kb*2*ndb + mt*kb*2*ndb + memcnt

!------------shape coefficient arrays and control volume metrics--------------------!

   ALLOCATE(a1u(0:nt,4))         ;a1u   = zero
   ALLOCATE(a2u(0:nt,4))         ;a2u   = zero 
   ALLOCATE(awx(0:nt,3))         ;awx   = zero 
   ALLOCATE(awy(0:nt,3))         ;awy   = zero 
   ALLOCATE(aw0(0:nt,3))         ;aw0   = zero 
   ALLOCATE(alpha(0:nt))         ;alpha = zero
   memcnt = nt*4*2*ndb + nt*3*3*ndb + nt*ndb + memcnt

!-----salinity and temperature bottom diffusion condition/bottom depth gradients----!

   ALLOCATE(phpn(0:mt))          ;phpn      = zero 
   ALLOCATE(pfpxb(mt))           ;pfpxb     = zero
   ALLOCATE(pfpyb(mt))           ;pfpyb     = zero
   ALLOCATE(sita_gd(0:mt))       ;sita_gd   = zero 
   ALLOCATE(ah_bottom(mt))       ;ah_bottom = zero 
   memcnt = mt*5*ndb + memcnt

   ALLOCATE(viscofh(0:mt,kb))    ;viscofh = zero
   ALLOCATE(viscofm(0:nt,kb))    ;viscofm = zero
   memcnt = mt*kb*ndb + nt*kb*ndb + memcnt

   ALLOCATE(hyw(0:mt,kb))        ;hyw     = zero
!-----special initialization which probably do nothing------------------------------!

   dt1(0)   = 100.0_sp

!---------------report approximate memory usage-------------------------------------!


   RETURN
   END SUBROUTINE alloc_vars
!==============================================================================|


END MODULE all_vars

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|

MODULE bcs
   USE mod_types
   USE mod_prec
   IMPLICIT NONE
   SAVE

!----------------boundary conditions: Julian tidal forcing--------------------------!

    TYPE(bc)                 :: elo_tm        !!TIME MAP FOR SURFACE ELEVATION DATA


!----------------GLobal Boundary Condition Information ----------------------------! 


   INTEGER,  ALLOCATABLE :: i_obc_gl(:)      !!GLOBAL ID OF OPEN BOUNDARY NODES
   INTEGER,  ALLOCATABLE :: i_obc_gl_w(:)    !!GLOBAL ID OF OPEN BOUNDARY NODES FOR WAVE
   INTEGER               :: iobcn_gl         !!LOCAL NUMBER OF OPEN BOUNDARY NODES FOR FVCOM
   INTEGER               :: iobcn_gl_w       !!LOCAL NUMBER OF OPEN BOUNDARY NODES FOR WAVE
   INTEGER               :: iobcn            !!LOCAL NUMBER OF OPEN BOUNDARY NODES FOR FVCOM
   INTEGER               :: iobcn_w          !!LOCAL NUMBER OF OPEN BOUNDARY NODES FOR SWAVE
   INTEGER,  ALLOCATABLE :: i_obc_n(:)       !!OPEN BOUNDARY NODE LIST FOR FVCOM
   INTEGER,  ALLOCATABLE :: i_obc_n_w(:)     !!OPEN BOUNDARY NODE LIST FOR SWAVE
   INTEGER,  ALLOCATABLE :: i_obc_n_output(:)!!LIST OF LOCAL OBC GLOBAL NODES FOR OUTPUT
   INTEGER,  ALLOCATABLE :: type_obc_gl(:)   !!OUTER BOUNDARY NODE TYPE (FOR SURFACE ELEVATION)
   INTEGER,  ALLOCATABLE :: type_obc(:)      !!OUTER BOUNDARY NODE TYPE (FOR SURFACE ELEVATION)
   INTEGER           ::  obc_ntime

!----------------boundary conditions: ground water----------------------------------!

  INTEGER, ALLOCATABLE :: node_bfw(:)      !!LOCAL GROUNDWATER NODES
  INTEGER, ALLOCATABLE :: bfw_gl2loc(:)    !!GLOBAL TO LOCAL MAPPING OF GWATER NODES
  REAL(sp),  ALLOCATABLE :: bfwqdis(:,:)   !!GROUNDWATER FRESH WATER FLUX DATA
!  TYPE(BC)      :: BFW_TM                  !!TIME MAP FOR GROUNDWATER DATA

!----------------boundary conditions: spectral tidal forcing----------------------!
   INTEGER :: ntidecomps
   REAL(sp),  ALLOCATABLE :: period(:)     !!TIDE PERIOD
   REAL(sp),  ALLOCATABLE :: apt(:,:)      !!TIDE AMPLITUDE
   REAL(sp),  ALLOCATABLE :: phai(:,:)     !!TIDE PHASE 
   REAL(sp),  ALLOCATABLE :: emean(:)      !!MEAN SURFACE ELEVATION

   ! ONLY USED IF COMPILED WITH EQUI_TIDE
   REAL(sp),  ALLOCATABLE :: apt_eqi(:)    !! EQUILIBIRUIM TIDE AMPLITUDE
    REAL(sp),  ALLOCATABLE :: beta_eqi(:)   !! EQUILIBIRUIM TIDE LOVE NUMBER
   CHARACTER(LEN=80), ALLOCATABLE :: tide_type(:)    !! EQUILIBIRUIM TIDE AMPLITUDE

   CHARACTER(LEN=12), PARAMETER :: diurnal="DIURNAL"
   CHARACTER(LEN=12), PARAMETER :: semidiurnal="SEMIDIURNAL"


!-- Old Tidal Periods before they became part of the forcing file ---------:    
!              s2         m2           n2          k1          p1         o1 
!PERIOD = (/43200.0_SP, 44712.0_SP, 45570.0_SP, 86164.0_SP, 86637.0_SP, 92950.0_SP/)


!----------------boundary conditions: Julian tidal forcing--------------------------!

   REAL(sp), ALLOCATABLE    :: elsbc(:,:)    !!INPUT SURFACE ELEVATION

END MODULE  bcs