mod_startup.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$
!/===========================================================================/

MODULE mod_startup
  USE mod_utils
  USE mod_nctools
  USE mod_input
  USE all_vars
  USE eqs_of_state
  USE mod_wd
  USE sinter



  


  
  IMPLICIT NONE
  
  PRIVATE
  PUBLIC :: read_ssh
  PUBLIC :: read_uv
  PUBLIC :: read_turb
  PUBLIC :: read_ts
  PUBLIC :: read_wetdry

  PUBLIC :: startup
!# if defined (DATA_ASSIM)
!  PUBLIC :: STARTUP_ASSIM
!# endif 
CONTAINS
!==============================================================================!
  SUBROUTINE startup
    IMPLICIT NONE

    IF(dbg_set(dbg_sbr)) write(ipt,*) "Start: Startup"

    IF(dbg_set(dbg_log)) THEN
       WRITE(ipt,*  )'!'
       WRITE(ipt,*  )'!           SETTING INITIAL CONDITIONS '
       WRITE(ipt,*  )'!'
    END IF

    IF (ASSOCIATED(nc_start))THEN
       IF(dbg_set(dbg_log)) write(ipt,*) "! NC_START FILE INFO:"
       CALL print_file(nc_start)
    END IF
       
!    SET THE SEA SURFACE HEIGHT
    if (dbg_set(dbg_log)) write(ipt,*) &
         & "! INITIALIZING SEA SURFACE HEIGHT"
    SELECT CASE(startup_type)
! =================================================
    CASE(startup_type_hotstart)
! =================================================

       CALL read_ssh

       IF(wetting_drying_on) CALL read_wetdry





! ------- New: Karsten Lettmann, June 2016 -----------------------
! ---------- end  new --------------------------------------------

!    if (dbg_set(dbg_log)) write(ipt,*) &
!         & "! INITIALIZING SEA ICE"



!    if (dbg_set(dbg_log)) write(ipt,*) &
!         & "! INITIALIZING SEA ICE222"





! =================================================
    CASE(startup_type_crashrestart)
! =================================================
       CALL read_ssh
     
       IF(wetting_drying_on) CALL read_wetdry





! ------- New: Karsten Lettmann, June 2016 -----------------------
! ---------- end  new --------------------------------------------






! =================================================
    CASE(startup_type_coldstart)
! =================================================

       CALL set_water_depth
       IF(wetting_drying_on) CALL set_wd_data

! =================================================
    CASE DEFAULT
! =================================================
       CALL fatal_error("STARTUP: UNKNOWN STARTUP TYPE: "//trim(startup_type))
    END SELECT


    if (dbg_set(dbg_log)) write(ipt,*) &
         & "! INITIALIZING VELOCITY FIELDS"

! SET STARTUP VALUES FOR VELOCITY
    SELECT CASE (startup_uv_type)
    CASE (startup_type_observed)
       CALL fatal_error("I DON'T KNOW HOW TO DO THAT KIND OF STARTUP")
    CASE(startup_type_linear) 
       CALL fatal_error("I DON'T KNOW HOW TO DO THAT KIND OF STARTUP")
    CASE(startup_type_constant)
       !CALL FATAL_ERROR("I DON'T KNOW HOW TO DO THAT KIND OF STARTUP")
       CALL set_constant_uv
    CASE(startup_type_default)
       ! OKAY - it is zero
    CASE(startup_type_setvalues)

       CALL read_uv

    CASE DEFAULT
       CALL fatal_error("UNKNOWN STARTUP_UV_TYPE")
    END SELECT

    if (dbg_set(dbg_log)) write(ipt,*) &
         & "! INITIALIZING TURBULENCE FIELDS"

! SET STARTUP VALUES FOR TURBULENCE
    SELECT CASE(startup_turb_type)
    CASE(startup_type_observed) 
       CALL fatal_error("I DON'T KNOW HOW TO DO THAT KIND OF STARTUP")
    CASE(startup_type_linear)
       CALL fatal_error("I DON'T KNOW HOW TO DO THAT KIND OF STARTUP")
    CASE(startup_type_constant)
       CALL fatal_error("I DON'T KNOW HOW TO DO THAT KIND OF STARTUP")
    CASE(startup_type_default)

       CALL set_default_turb

    CASE(startup_type_setvalues)

       CALL read_turb

    CASE DEFAULT
       CALL fatal_error("UNKNOWN STARTUP_TURB_TYPE")
    END SELECT

    if (dbg_set(dbg_log)) write(ipt,*) &
         & "! INITIALIZING TEMPERATURE AND SALINITY"

! SET STARTUP VALUES FOR TEMPERATER AND SALINITY
    SELECT CASE(startup_ts_type)
    CASE(startup_type_observed) 

       CALL set_observed_ts

    CASE(startup_type_linear)

       CALL set_linear_ts

    CASE(startup_type_constant)

       CALL set_constant_ts

    CASE(startup_type_default)
       CALL fatal_error("There is no default startup for Temperature and Salinity")
    CASE(startup_type_setvalues)

       CALL read_ts

    CASE DEFAULT
       CALL fatal_error("UNKNOWN STARTUP_TS_TYPE")
    END SELECT



    ! ONCE WE HAVE STARTED THE MODEL POINT THE START FILE OBJECT TO
    ! ITS OWN OUTPUT TO RELOAD OLD TIME STATES
    IF(.not. associated(nc_start,nc_rst)) THEN
       
       ! SOME START UPS DO NOT HAVE A START FILE
       IF(Associated(nc_start)) CALL kill_file(nc_start)

       nc_start => nc_rst
    END IF

    IF(dbg_set(dbg_sbr)) write(ipt,*) "End: Startup"


  END SUBROUTINE startup
!==============================================================================!

!# if defined (DATA_ASSIM)
!!==============================================================================!
!  SUBROUTINE STARTUP_ASSIM
!# if defined(SEDIMENT) 
!  USE MOD_SED, ONLY: sed_hot_start
!# endif
!    IMPLICIT NONE

!    IF(DBG_SET(DBG_SBR)) write(ipt,*) "Start: Startup_ASSIM"
       
!!    SET THE SEA SURFACE HEIGHT
!    if (dbg_set(dbg_log)) write(ipt,*) &
!         & "! INITIALIZING SEA SURFACE HEIGHT"
!       CALL READ_SSH

!       IF(WETTING_DRYING_ON) CALL READ_WETDRY

!# if defined (NH)
!       CALL READ_NH
!# endif

!# if defined (ATMO_TIDE)
!       CALL READ_ATMO
!# endif

!# if defined (EQUI_TIDE)
!       CALL READ_EQI
!# endif

!# if defined (ICE)    
!       CALL  READ_ICE
!       CALL  AGGREGATE  !!  ggao 07312008
!# endif    

!# if defined (SEDIMENT)
!  if(sed_hot_start)then
!    CALL READ_SED
!  endif
!# endif


!       CALL READ_UV

!       CALL READ_TURB

!       CALL READ_TS
!    if (dbg_set(dbg_log)) write(ipt,*) &
!         & "! INITIALIZING TS"


!#   if defined (WATER_QUALITY)

!       CALL READ_WQM

!#   endif

!#   if defined (BioGen)

!       CALL READ_BIO

!#   endif


!    IF(DBG_SET(DBG_SBR)) write(ipt,*) "End: Startup_Assim"



!  END SUBROUTINE STARTUP_ASSIM
!!==============================================================================!
!# endif

  SUBROUTINE read_ssh
    IMPLICIT NONE
    TYPE(ncvar),  POINTER :: var
    TYPE(ncdim),  POINTER :: dim
    LOGICAL :: found
    INTEGER :: stkcnt

    IF(dbg_set(dbg_sbr)) write(ipt,*) "Start: READ_SSH"
    
    stkcnt = nc_start%FTIME%PREV_STKCNT

    ! LOAD EL
    var => find_var(nc_start,'zeta',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'zeta'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, el)
    CALL nc_read_var(var,stkcnt)

    ! LOAD ET
    var => find_var(nc_start,'et',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'et'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, et)
    CALL nc_read_var(var,stkcnt)

    !----------------------------------------------------------------
    ! Read the most recent bathymetry if Morphodynamics is Active
    !----------------------------------------------------------------


    !----------------------------------------------------------------
    ! Given SSH and Bathy, Update the Bathymetry 
    !----------------------------------------------------------------
    d  = h + el
    dt = h + et

    
    CALL n2e2d(h,h1)
    CALL n2e2d(el,el1)
    CALL n2e2d(d,d1)
    CALL n2e2d(dt,dt1)

    IF(dbg_set(dbg_sbr)) write(ipt,*) "Start: READ_SSH"

  END SUBROUTINE read_ssh
!==============================================================================!
!==============================================================================!
  SUBROUTINE read_wetdry
    USE mod_wd
    IMPLICIT NONE
    TYPE(ncvar),  POINTER :: var
    TYPE(ncdim),  POINTER :: dim
    LOGICAL :: found
    INTEGER :: stkcnt

    IF(dbg_set(dbg_sbr)) write(ipt,*) "Start: READ_WETDRY"

    stkcnt = nc_start%FTIME%PREV_STKCNT

    ! LOAD ISWETN
    var => find_var(nc_start,'wet_nodes',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'wet_nodes'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, iswetn)
    CALL nc_read_var(var,stkcnt)

    ! LOAD ISWETC
    var => find_var(nc_start,'wet_cells',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'wet_cells'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, iswetc)
    CALL nc_read_var(var,stkcnt)


    ! LOAD ISWETNT
    var => find_var(nc_start,'wet_nodes_prev_int',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'wet_nodes_prev_int'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, iswetnt)
    CALL nc_read_var(var,stkcnt)

    ! LOAD ISWETCT
    var => find_var(nc_start,'wet_cells_prev_int',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'wet_cells_prev_int'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, iswetct)
    CALL nc_read_var(var,stkcnt)

    ! LOAD ISWETCE
    var => find_var(nc_start,'wet_cells_prev_ext',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'wet_cells_prev_ext'&
         & IN THE HOTSTART FILE OBJECT")
    ! ------ new: Karsten Lettmann, May 2016 -------------
    !CALL NC_CONNECT_AVAR(VAR, ISWETC)  ! original line
    CALL nc_connect_avar(var, iswetce)
    ! ------------- end new ------------------------------
    CALL nc_read_var(var,stkcnt)

    IF(dbg_set(dbg_sbr)) write(ipt,*) "End: READ_WETDRY"
    
  END SUBROUTINE read_wetdry
!==============================================================================!
!==============================================================================!
  SUBROUTINE read_atmo
    IMPLICIT NONE
    TYPE(ncvar),  POINTER :: VAR
    TYPE(ncdim),  POINTER :: DIM
    LOGICAL :: FOUND
    INTEGER :: STKCNT
    
    IF(dbg_set(dbg_sbr)) write(ipt,*) "START: READ_ATMO" 

    stkcnt = nc_start%FTIME%PREV_STKCNT

    ! LOAD ATMO
    var => find_var(nc_start,'el_atmo',found)
    IF(.not. found) CALL fatal_error&
         &("COULD NOT FIND VARIABLE 'el_atmo' IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, el_atmo)
    CALL nc_read_var(var,stkcnt)

    IF(dbg_set(dbg_sbr)) write(ipt,*) "End: READ_ATMO"

  END SUBROUTINE read_atmo
!==============================================================================!
  SUBROUTINE read_eqi
    IMPLICIT NONE
    TYPE(ncvar),  POINTER :: VAR
    TYPE(ncdim),  POINTER :: DIM
    LOGICAL :: FOUND
    INTEGER :: STKCNT

    IF(dbg_set(dbg_sbr)) write(ipt,*) "Start: READ_EQI" 
    
    stkcnt = nc_start%FTIME%PREV_STKCNT

    ! LOAD ATMO
    var => find_var(nc_start,'el_eqi',found)
    IF(.not. found) CALL fatal_error&
         &("COULD NOT FIND VARIABLE 'el_eqi' IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, el_eqi)
    CALL nc_read_var(var,stkcnt)

    IF(dbg_set(dbg_sbr)) write(ipt,*) "End: READ_EQI" 
   
  END SUBROUTINE read_eqi
!==============================================================================!

! ------- New: Karsten Lettmann, June 2016 -----------------------
! ------------- end new ---------------------------------------------------------

!==============================================================================!
  SUBROUTINE read_ts
    IMPLICIT NONE
    TYPE(ncvar),  POINTER :: var
    TYPE(ncdim),  POINTER :: dim
    LOGICAL :: found
    INTEGER :: stkcnt, k
    REAL(sp), DIMENSION(0:MT,KB) :: pressure
    
    IF(dbg_set(dbg_sbr)) write(ipt,*) "Start: READ_TS" 
   
    stkcnt = nc_start%FTIME%PREV_STKCNT


    ! LOAD TEMPERATURE
    var => find_var(nc_start,'temp',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'temp'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, t1)
    CALL nc_read_var(var,stkcnt)


    ! LOAD MEAN INITIAL TEMPERATURE
    var => find_var(nc_start,'tmean1',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'tmean1'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, tmean1)
    CALL nc_read_var(var)

    ! LOAD SALINITY
    var => find_var(nc_start,'salinity',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'saltinity'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, s1)
    CALL nc_read_var(var,stkcnt)

    ! LOAD MEAN INITIAL SALINITY
    var => find_var(nc_start,'smean1',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'smean1'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, smean1)
    CALL nc_read_var(var)

    ! LOAD DENSITY
    var => find_var(nc_start,'rho1',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'rho1'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, rho1)
    CALL nc_read_var(var,stkcnt)

    ! LOAD MEAN DENSITY
    var => find_var(nc_start,'rmean1',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'rmean1'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, rmean1)
    CALL nc_read_var(var)


    ! AVERAGE FROM CELLS TO FACE CENTERS


!JQI    SELECT CASE(SEA_WATER_DENSITY_FUNCTION)
!JQI    CASE(SW_DENS1)

!JQI       ! SET MEAN DENSITY
!JQI       DO K=1,KBM1
!JQI          PRESSURE(:,K) = -GRAV_N*1.025_SP*(ZZ(:,K)*D(:))*0.1_SP
!JQI       END DO
!JQI       CALL FOFONOFF_MILLARD(SMEAN1,TMEAN1,PRESSURE,0.0_SP,RMEAN1)
!JQI       RMEAN1(0,:)=0.0_SP
!JQI       RMEAN1(:,KB)=0.0_SP

!JQI       ! SET REAL DENSITY
!JQI       CALL DENS1 ! GENERIC CALL TO FOFONOFF_MILLARD FOR S1,T1...
       
!JQI    CASE(SW_DENS2)
!JQI       ! SET MEAN DENSITY
!JQI       CALL DENS2G(SMEAN1,TMEAN1,RMEAN1)
!JQI       RMEAN1(0,:)=0.0_SP
!JQI       RMEAN1(:,KB)=0.0_SP

!JQI       ! SET REAL DENSITY
!JQI       CALL DENS2 ! GENERIC CALL TO DENS2G FOR S1,T1...

!JQI    CASE(SW_DENS3)

!JQI       ! SET MEAN DENSITY
!JQI       DO K=1,KBM1
!JQI          PRESSURE(:,K) = -GRAV_N*1.025_SP*(ZZ(:,K)*D(:))*0.01_SP
!JQI       END DO
!JQI       CALL JACKET_MCDOUGALL(SMEAN1,TMEAN1,PRESSURE,RMEAN1)
!JQI       RMEAN1(0,:)=0.0_SP
!JQI       RMEAN1(:,KB)=0.0_SP

!JQI       ! SET REAL DENSITY
!JQI       CALL DENS3 ! GENERIC CALL TO JACKET_MCDOUGALL FOR S1,T1.

!JQI    CASE DEFAULT
!JQI       CALL FATAL_ERROR("INVALID DENSITY FUNCTION SELECTED:",&
!JQI            & "   "//TRIM(SEA_WATER_DENSITY_FUNCTION) )
!JQI    END SELECT

!J. GE
    t0=t1
    s0=s1
!J. GE    
    CALL n2e3d(t1,t)
    CALL n2e3d(s1,s)
    CALL n2e3d(rho1,rho)
    CALL n2e3d(tmean1,tmean)
    CALL n2e3d(smean1,smean)
    CALL n2e3d(rmean1,rmean)

   
   IF(dbg_set(dbg_sbr)) write(ipt,*) "End: READ_TS" 
   
 END SUBROUTINE read_ts
!==============================================================================!
  SUBROUTINE set_observed_ts
    IMPLICIT NONE
    TYPE(ncvar),  POINTER :: var
    TYPE(ncdim),  POINTER :: dim
    LOGICAL :: found
    INTEGER :: stkcnt
    INTEGER ksl                                 !!NUMBER OF STANDARD SEA LEVELS 
    REAL(sp), ALLOCATABLE, TARGET :: dpthsl(:)          !!DEPTH AT STANDARD SEA LEVEL
    REAL(sp), ALLOCATABLE, TARGET :: tsl(:,:),ssl(:,:)  !!T/S AT STANDARD SEA LEVEL
    REAL(sp),ALLOCATABLE          :: ta(:),sa(:)

    REAL(sp),DIMENSION(KBM1)      :: ti,si,zi
    REAL(sp),DIMENSION(0:MT,KB) :: pressure
    INTEGER :: k, i, ib, ierr
    REAL(sp) :: scount, fac, rbuf
    
    IF(dbg_set(dbg_sbr)) write(ipt,*) "START: SET_OBSERVED_TS" 

    ti       = 0.0_sp
    si       = 0.0_sp
    zi       = 0.0_sp
    pressure = 0.0_sp

    stkcnt = nc_start%FTIME%PREV_STKCNT

    dim => find_dim(nc_start,'ksl',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND DIMENSION 'ksl'&
         & IN THE STARTUP FILE OBJECT")

    ksl = dim%DIM
    ! ALLOCATE SPACE
    ALLOCATE(tsl(mt,ksl))
    ALLOCATE(ssl(mt,ksl))
    ALLOCATE(dpthsl(ksl))
    tsl   = 0.0_sp
    ssl   = 0.0_sp
    dpthsl= 0.0_sp

    ! ALLOCATE SPACE FOR THE AVERAGE TEMPERATURE AT EACH DEPTH
    ALLOCATE(ta(ksl),sa(ksl))
    ta = 0.0_sp
    sa = 0.0_sp

    ! READ THE DATA
    var => find_var(nc_start,'tsl',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'tsl'&
         & IN THE STARTUP FILE OBJECT")
    CALL nc_connect_avar(var, tsl)
    CALL nc_read_var(var,stkcnt)
    CALL nc_disconnect(var)

    var => find_var(nc_start,'ssl',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'ssl'&
         & IN THE STARTUP FILE OBJECT")
    CALL nc_connect_avar(var, ssl)
    CALL nc_read_var(var,stkcnt)
    CALL nc_disconnect(var)


    var => find_var(nc_start,'zsl',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'zsl'&
         & IN THE STARTUP FILE OBJECT")
    CALL nc_connect_avar(var, dpthsl)
    CALL nc_read_var(var)
    CALL nc_disconnect(var)


!==============================================================================|
!    CALCULATE HORIZONTAL AVERAGE VALUES OF TEMPERATURE/SALINITY/DENSITY       !
!==============================================================================|
   
      where(ssl <0.0_sp) ssl=0.0
 
    ! GET THE AVERAGE T/S AT EACH DEPTH ON STANDARD LEVELS
    IF(serial)THEN
       DO k = 1, ksl
          scount = 0.0_sp
          DO i = 1, mt
             IF(-h(i) <= dpthsl(k)) THEN
                scount = scount + 1.0_sp
                ta(k) = ta(k) + tsl(i,k)
                sa(k) = sa(k) + ssl(i,k)
             END IF
          END DO
          IF(scount >= 1.0_sp)THEN
             ta(k)=ta(k)/scount
             sa(k)=sa(k)/scount
          END IF
       END DO
    END IF
    
    
    
    ! NOW INTERPOLATE FROM STANDARD LEVELS TO THE VALUE AT EACH NODE
    DO i=1,mt
       DO k=1,kbm1
          zi(k)=zz(i,k)*d(i)+el(i)
       END DO
       
       ! LEVEL AVERAGE T AND S
       CALL sinter_extrp_up(dpthsl,ta,zi,ti,ksl,kbm1)
       CALL sinter_extrp_up(dpthsl,sa,zi,si,ksl,kbm1)
       
       tmean1(i,1:kbm1) = ti(:)
       smean1(i,1:kbm1) = si(:)
              
       ! REAL T AND S
       CALL sinter_extrp_up(dpthsl,tsl(i,:),zi,ti,ksl,kbm1)
       CALL sinter_extrp_up(dpthsl,ssl(i,:),zi,si,ksl,kbm1)
       
       t1(i,1:kbm1) = ti(:)
       s1(i,1:kbm1) = si(:)
!J. Ge
       t0(i,1:kbm1) = ti(:)
       s0(i,1:kbm1) = si(:)
!J. Ge
    END DO

       where(s1<0.0_sp) s1=0.0

    IF(.NOT.barotropic)THEN
     SELECT CASE(sea_water_density_function)
     CASE(sw_dens1)

       ! SET MEAN DENSITY
       DO k=1,kbm1
          pressure(:,k) = -grav_n*1.025_sp*(zz(:,k)*d(:))*0.1_sp
       END DO
       CALL fofonoff_millard(smean1,tmean1,pressure,0.0_sp,rmean1)
       rmean1(0,:)=0.0_sp
       rmean1(:,kb)=0.0_sp

       ! SET REAL DENSITY
       CALL dens1 ! GENERIC CALL TO FOFONOFF_MILLARD FOR S1,T1...

     CASE(sw_dens2)

       ! SET MEAN DENSITY
       CALL dens2g(smean1,tmean1,rmean1)
       rmean1(0,:)=0.0_sp
       rmean1(:,kb)=0.0_sp

       ! SET REAL DENSITY
       CALL dens2 ! GENERIC CALL TO DENS2G FOR S1,T1...

     CASE(sw_dens3)
     
      ! SET MEAN DENSITY
       DO k=1,kbm1
          pressure(:,k) = -grav_n*1.025_sp*(zz(:,k)*d(:))*0.01_sp
       END DO
       CALL jacket_mcdougall(smean1,tmean1,pressure,rmean1)
       rmean1(0,:)=0.0_sp
       rmean1(:,kb)=0.0_sp

       ! SET REAL DENSITY
       CALL dens3 ! GENERIC CALL TO JACKET_MCDOUGALL FOR S1,T1.

     CASE DEFAULT
      CALL fatal_error("INVALID DENSITY FUNCTION SELECTED:",&
           & "   "//trim(sea_water_density_function) )
     END SELECT
    ELSE
     rho1   = 2.3e-2_sp
     rho    = 2.3e-2_sp
     rmean1 = 2.3e-2_sp
    END IF 
   

   CALL n2e3d(t1,t)
   CALL n2e3d(s1,s)
   CALL n2e3d(tmean1,tmean)
   CALL n2e3d(smean1,smean)
   CALL n2e3d(rmean1,rmean)
   ! DENSITY IS ALREADY INTERPOLATED TO ELEMENTS FOR RHO IN DENSX

   DEALLOCATE(tsl,ssl,dpthsl)
   DEALLOCATE(ta,sa)


   IF(dbg_set(dbg_sbr)) write(ipt,*) "END: SET_OBSERVED_TS" 
    
  END SUBROUTINE set_observed_ts
!==============================================================================!
  SUBROUTINE set_linear_ts
    IMPLICIT NONE

    ! BY DEFINITION OF LINEAR...
    INTEGER, PARAMETER :: ksl =2

    REAL(sp), DIMENSION(KBM1)   :: ti,si,zi
    REAL(sp), DIMENSION(KSL)   :: ta,sa,za
    INTEGER :: k, i

    IF(dbg_set(dbg_sbr)) write(ipt,*) "START: SET_LINEAR_TS" 

    za(1) = 0.0_sp
    za(2) = startup_dmax
    
    ta(1) = startup_t_vals(1)
    ta(2) = startup_t_vals(2)
    
    sa(1) = startup_s_vals(1)
    sa(2) = startup_s_vals(2)
    
    ! THE HORIZONTAL AVERAGE VALUE IS ALSO THE TRUE VALUE SINCE THE
    ! LINEAR EQUATION DOES NOT VARY WITH LOCATION 

    DO i = 1, mt
       DO k=1,kbm1
          ! CALCULATE ZI relative to z=0
          zi(k)=zz(i,k)*d(i)+el(i)
       END DO
    
       CALL sinter_extrp_up(za,ta,zi,ti,ksl,kbm1)
       CALL sinter_extrp_up(za,sa,zi,si,ksl,kbm1)
    
       t1(i,1:kbm1) = ti
       s1(i,1:kbm1) = si
!J. Ge
       t0(i,1:kbm1) = ti
       s0(i,1:kbm1) = si
!J. Ge

    END DO

    IF(.NOT.barotropic)THEN
     SELECT CASE(sea_water_density_function)
     CASE(sw_dens1)
       CALL dens1 ! USE GENERIC INTERFACE TO DENSX
     CASE(sw_dens2)
       CALL dens2
     CASE(sw_dens3)
       CALL dens3
     CASE DEFAULT
       CALL fatal_error("INVALID DENSITY FUNCTION SELECTED:",&
           & "   "//trim(sea_water_density_function) )
     END SELECT
    ELSE
     rho1 = 2.3e-2_sp
     rho  = 2.3e-2_sp
    END IF 
   
   ! FOR LINEAR AND CONSTANT, THE MEAN IS EQUAL TO THE INITIAL VALUE

   tmean1=t1
   smean1=s1
   rmean1=rho1
   rmean=rho
   CALL n2e3d(t1,t)
   CALL n2e3d(s1,s)
   tmean=t
   smean=s
   
    IF(dbg_set(dbg_sbr)) write(ipt,*) "END: SET_LINEAR_TS" 
    
  END SUBROUTINE set_linear_ts
!==============================================================================!
  SUBROUTINE set_constant_ts
    IMPLICIT NONE
    
    IF(dbg_set(dbg_sbr)) write(ipt,*) "START: SET_CONSTANT_TS" 


    t1(:,1:kbm1) = startup_t_vals(1)
    s1(:,1:kbm1) = startup_s_vals(1)
!J. Ge
    t0(:,1:kbm1) = startup_t_vals(1)
    s0(:,1:kbm1) = startup_s_vals(1)
!J. Ge

    IF(.NOT.barotropic)THEN
     SELECT CASE(sea_water_density_function)
     CASE(sw_dens1)
       CALL dens1
     CASE(sw_dens2)
       CALL dens2
     CASE(sw_dens3)
       CALL dens3
     CASE DEFAULT
       CALL fatal_error("INVALID DENSITY FUNCTION SELECTED:",&
            & "   "//trim(sea_water_density_function) )
     END SELECT
    ELSE
     rho1 = 2.3e-2_sp
     rho  = 2.3e-2_sp
    END IF 

   ! FOR LINEAR AND CONSTANT, THE MEAN IS EQUAL TO THE INITIAL VALUE

   tmean1=t1
   smean1=s1
   rmean1=rho1

   CALL n2e3d(t1,t)
   CALL n2e3d(s1,s)
   
   tmean=t
   rmean=rho
   smean=s

   IF(dbg_set(dbg_sbr)) write(ipt,*) "END: SET_CONSTANT_TS" 

 END SUBROUTINE set_constant_ts
!==============================================================================!
!==============================================================================!
  SUBROUTINE set_constant_uv
    IMPLICIT NONE

    IF(dbg_set(dbg_sbr)) write(ipt,*) "START: SET_CONSTANT_UV"


    u(:,1:kbm1) = startup_u_vals    
    v(:,1:kbm1) = startup_v_vals
    ua          = startup_u_vals
    va          = startup_v_vals
    wts         = 0.0          
    w           = 0.0

   IF(dbg_set(dbg_sbr)) write(ipt,*) "END: SET_CONSTANT_UV"

 END SUBROUTINE set_constant_uv
!==============================================================================!

  SUBROUTINE read_uv
    IMPLICIT NONE
    TYPE(ncvar),  POINTER :: var
    TYPE(ncdim),  POINTER :: dim
    LOGICAL :: found
    INTEGER :: stkcnt

    IF(dbg_set(dbg_sbr)) write(ipt,*) "START: READ_UV" 


    stkcnt = nc_start%FTIME%PREV_STKCNT

    var => find_var(nc_start,'u',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'u'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, u)
    CALL nc_read_var(var,stkcnt)


    var => find_var(nc_start,'v',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'v'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, v)
    CALL nc_read_var(var,stkcnt)

    var => find_var(nc_start,'omega',found)
    IF(found) THEN
       CALL nc_connect_avar(var, wts)
       CALL nc_read_var(var,stkcnt)
       
       CALL n2e3d(wts,w)
    ELSE
       var => find_var(nc_start,'w',found)
       IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'w' &
            & or 'omega' IN THE HOTSTART FILE OBJECT")
       CALL nc_connect_avar(var, w)
       CALL nc_read_var(var,stkcnt)
    END IF

    var => find_var(nc_start,'ua',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'ua'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, ua)
    CALL nc_read_var(var,stkcnt)

    var => find_var(nc_start,'va',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'va'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, va)
    CALL nc_read_var(var,stkcnt)

    IF(dbg_set(dbg_sbr)) write(ipt,*) "END: READ_UV" 

  END SUBROUTINE read_uv
!==============================================================================!
  SUBROUTINE read_turb
    IMPLICIT NONE
    TYPE(ncvar),  POINTER :: var
    TYPE(ncdim),  POINTER :: dim
    LOGICAL :: found
    INTEGER :: stkcnt

    IF(dbg_set(dbg_sbr)) write(ipt,*) "START: READ_TURB" 

    stkcnt = nc_start%FTIME%PREV_STKCNT


    var => find_var(nc_start,'q2',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'q2'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, q2)
    CALL nc_read_var(var,stkcnt)

    var => find_var(nc_start,'q2l',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'q2l'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, q2l)
    CALL nc_read_var(var,stkcnt)

    var => find_var(nc_start,'l',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'l'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, l)
    CALL nc_read_var(var,stkcnt)


    var => find_var(nc_start,'km',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'km'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, km)
    CALL nc_read_var(var,stkcnt)

    var => find_var(nc_start,'kq',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'kq'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, kq)
    CALL nc_read_var(var,stkcnt)

    var => find_var(nc_start,'kh',found)
    IF(.not. found) CALL fatal_error("COULD NOT FIND VARIABLE 'kh'&
         & IN THE HOTSTART FILE OBJECT")
    CALL nc_connect_avar(var, kh)
    CALL nc_read_var(var,stkcnt)


    CALL n2e3d(km,km1)

    IF(dbg_set(dbg_sbr)) write(ipt,*) "END: READ_TURB" 

  END SUBROUTINE read_turb
!==============================================================================|
   SUBROUTINE set_default_turb         
!==============================================================================|
!   Initialize Turbulent Kinetic Energy and Length Scale                       |
!==============================================================================|
   IMPLICIT NONE
   INTEGER :: i,k
!==============================================================================|
   IF(dbg_set(dbg_sbr)) write(ipt,*) "START: SET_DEFAULT_TURB" 
   
   DO k=1,kbm1
        aam(:,k) = nn_hvc(:)
   END DO

!
!------------------------BOUNDARY VALUES---------------------------------------!
!

   DO i = 1, mt
     km(i,1)   = 0.0_sp
     km(i,kb)  = 0.0_sp
     kh(i,1)   = 0.0_sp
     kh(i,kb)  = 0.0_sp
     kq(i,1)   = 0.0_sp
     kq(i,kb)  = 0.0_sp
     l(i,1)    = 0.0_sp
     l(i,kb)   = 0.0_sp
     q2(i,1)   = 0.0_sp
     q2(i,kb)  = 0.0_sp
     q2l(i,1)  = 0.0_sp
     q2l(i,kb) = 0.0_sp
   END DO


!
!------------------------INTERNAL VALUES---------------------------------------!
!
   DO  k = 2, kbm1
     DO i = 1, mt
       IF (d(i) > 0.0_sp) THEN
         q2(i,k)  = 1.e-8
         q2l(i,k) = 1.e-8
         l(i,k)   = 1.
         km(i,k)  = 2.*umol
         kq(i,k)  = 2.*umol
         kh(i,k)  = 2.*umol / vprnu
       END IF
     END DO
   END DO
   CALL n2e3d(km,km1)

   IF(dbg_set(dbg_sbr)) write(ipt,*) "END: SET_DEFAULT_TURB" 

 END SUBROUTINE set_default_turb
!==============================================================================|

!==============================================================================|
!   READ IN STATIC WATER DEPTH AND CALCULATE RELATED QUANTITIES                |
!                                                                              |
!   INPUTS: H(NNODE) BATHYMETRIC DEPTH AT NODES                    |
!   INITIALIZES: D(NNODE) DEPTH AT NODES                       |
!   INITIALIZES: DT(NNODE) ???                                 |
!   INITIALIZES: H1(NNODE) BATHYMETRIC DEPTH AT ELEMENTS               |
!   INITIALIZES: D1(NNODE) DEPTH AT NODES                              | 
!   INITIALIZES: DT1(NNODE) ??                                             |
!==============================================================================|

   SUBROUTINE set_water_depth         
!------------------------------------------------------------------------------|
   USE mod_obcs
   USE mod_nesting
   IMPLICIT NONE
   REAL(sp) :: temp
   INTEGER  :: i,k,j1,j2
!------------------------------------------------------------------------------|
   IF(dbg_set(dbg_sbr)) write(ipt,*) "START: SET_WATER_DEPTH" 


!
!  ADJUST STATIC HEIGHT AND CALCULATE DYNAMIC DEPTHS (D) AND (DT)
!
   h  = h + static_ssh_adj
   d  = h + el
   dt = h + et

!
!  ADJUST SIGMA VALUES ON OUTER BOUNDARY 
!

   IF(iobcn > 0 .AND. obc_depth_control_on) THEN
     IF(.NOT. nesting_on)THEN
     DO i = 1,iobcn
       j1 = i_obc_n(i)
       j2 = next_obc(i)
       h(j1) = h(j2)
       d(j1) = d(j2)
       dt(j1) = dt(j2)
       
       DO k = 1,kb
         z(j1,k)   = z(j2,k)
         zz(j1,k)  = zz(j2,k)
         dz(j1,k)  = dz(j2,k)
         dzz(j1,k) = dzz(j2,k)
       END DO    
     END DO
     END IF
   END IF

! 
!  CALCULATE FACE-CENTERED VALUES OF BATHYMETRY AND DEPTH
!
   DO i=1,nt    
     h1(i)  = (h(nv(i,1))+h(nv(i,2))+h(nv(i,3)))/3.0_sp
     d1(i)  = h1(i)+el1(i)
     dt1(i) = h1(i)+et1(i)
     
     DO k = 1,kb
       z1(i,k)=(z(nv(i,1),k)+z(nv(i,2),k)+z(nv(i,3),k))/3.0_sp
     END DO  
   END DO

!-----AFTER MODIFYING BOUNDARY SIGMA VALUES ------------------------------------!
!-----RECOMPUTE SIGMA DERIVATIVES AND INTRA SIGMA LEVELS AGAIN ON CELL----------!

   DO k=1,kb-1
     DO i=1,nt
       dz1(i,k)  = z1(i,k)-z1(i,k+1)
       zz1(i,k)  = .5_sp*(z1(i,k)+z1(i,k+1))
     END DO
   END DO

   DO i=1,nt
     zz1(i,kb) = 2.0_sp*zz1(i,kb-1)-zz1(i,kb-2)
   END DO

   DO k=1,kbm2
     DO i=1,nt
       dzz1(i,k) = zz1(i,k)-zz1(i,k+1)
     END DO
   END DO
   
   dzz1(:,kbm1) = 0.0_sp
   dz1(:,kb)    = 0.0_sp

   IF(dbg_set(dbg_sbr)) write(ipt,*) "END: SET_WATER_DEPTH" 
   RETURN
 END SUBROUTINE set_water_depth
!==============================================================================|








END MODULE mod_startup