mod_setup.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_setup
  USE all_vars
  USE mod_par
  USE mod_input
  USE mod_ncdio
  USE mod_nctools
  USE mod_utils
  USE mod_obcs
  USE mod_force

  IMPLICIT NONE

  SAVE
  ! COORDINATE VARIABLES
  REAL(sp), ALLOCATABLE, TARGET :: x_gbl(:),y_gbl(:)
  REAL(sp), ALLOCATABLE, TARGET :: x_lcl(:),y_lcl(:)
  
  ! CORIOLIS VARIABLES
  REAL(sp), ALLOCATABLE, TARGET :: c_lcl(:)     !!CORIOLIS PARAMETER AT NODES

  ! DEPTH VARIABLE
  REAL(sp), ALLOCATABLE, TARGET :: h_lcl(:)     !! DEPTH PARAMETER AT NODES

  ! SPONGE LAYER VARIABLE
  INTEGER, ALLOCATABLE, TARGET  :: n_spg(:)
  REAL(sp), ALLOCATABLE, TARGET :: r_spg(:),c_spg(:),x_spg(:),y_spg(:)
  INTEGER :: nsponge

  PRIVATE :: sigma_geometric, sigma_generalized, sigma_tanh

CONTAINS

!==============================================================================!
  SUBROUTINE setup_center_coords
    USE mod_spherical
    IMPLICIT NONE
    INTEGER I,J,IERR,STATUS,SENDID
    REAL(SP) SBUF

    INTEGER K,ITMP
    REAL(DP) VX1,VY1,VX2,VY2,VX3,VY3
    REAL(DP) EVX12,EVX13,EVX23,EVY12,EVY13,EVY23,EVXY
    REAL(DP) VX12,VY12,VX23,VY23,VX31,VY31
    INTEGER :: SENDER

    ! David Added:
    REAL(SP), allocatable :: xc_buf(:), lonc_buf(:)
    REAL(SP), allocatable :: yc_buf(:), latc_buf(:)
 
    IF(dbg_set(dbg_sbr)) write(ipt,*) "SETUP_CENTER_COORDS: START"
    
    
!==============================================================================|
!   SET UP LOCAL MESH (HORIZONTAL COORDINATES)                                 |
!==============================================================================|
!--------------CALCULATE GLOBAL MINIMUMS AND MAXIMUMS--------------------------!
    
! IF THIS IS NOT A SPHERICAL MODEL REMOVE MIN COORDINATE VALUE    
    
    vxmin = minval(vx(1:mt)) ; vxmax = maxval(vx(1:mt))
    vymin = minval(vy(1:mt)) ; vymax = maxval(vy(1:mt))
    
    !--------------SHIFT GRID TO UPPER RIGHT CARTESIAN-----------------------------!
    
    vx = vx - vxmin
    vy = vy - vymin
    vx(0) = 0.0_sp ; vy(0) = 0.0_sp
    
    !--------------CALCULATE GLOBAL ELEMENT CENTER GRID COORDINATES----------------!
    CALL n2e2d(vx,xc)
    CALL n2e2d(vy,yc)
    
    xmc = xc + vxmin
    xmc(0)= 0.0_sp
    ymc = yc + vymin
    ymc(0)= 0.0_sp

    IF (use_proj ) THEN
        IF (serial) THEN
           CALL meters2degrees(xmc(1:nt),ymc(1:nt)&
                &   ,projection_reference,lonc(1:nt),latc(1:nt),nt)
        END IF

    END IF

    

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

  END SUBROUTINE setup_center_coords
  !==============================================================================|
  !   SET HORIZONTAL MIXING_COEFFICIENT                                      |
  !==============================================================================|
  SUBROUTINE setup_horizontal_mixing_coefficient
    IMPLICIT NONE
  
    if (horizontal_mixing_kind .eq. sttc) THEN

       if(dbg_set(dbg_log)) then
          write(ipt,*) "! "
          write(ipt,*) "! Setting Staticly Variable viscosity"
       end if

       CALL load_horizontal_mixing_coefficient(nn_hvc,cc_hvc)
       
       
    else if(horizontal_mixing_kind .eq. cnstnt) THEN
       
       cc_hvc=horizontal_mixing_coefficient
       cc_hvc(0)=0.0_sp
       nn_hvc=horizontal_mixing_coefficient
       nn_hvc(0)=0.0_sp
       
    else
       CALL fatal_error&
            &("HORIZONTAL_MIXING_KIND ERROR",&
            & "This should not happen")
       
    end if
    
    
  END SUBROUTINE setup_horizontal_mixing_coefficient
  !==============================================================================|
  !   SET BOTTOM ROUGHNESS                                                       |
  !==============================================================================|
  SUBROUTINE setup_bottom_roughness
    IMPLICIT NONE
  
    if (bottom_roughness_kind .eq. sttc) THEN
  
       if(dbg_set(dbg_log)) then
          write(ipt,*) "! "
          write(ipt,*) "! Setting Staticly Variable Bottom Roughness"
       end if

       CALL load_bottom_roughness(cc_z0b)
              
    else if(bottom_roughness_kind .eq. cnstnt) THEN
       
       cc_z0b    = bottom_roughness_lengthscale
       cc_z0b(0) = 0.0_sp
       
    else
       CALL fatal_error&
            &("HORIZONTAL_MIXING_KIND ERROR",&
            & "This should not happen")
       
    end if
    
    
  END SUBROUTINE setup_bottom_roughness
  !==============================================================================|
  !   SET UP LOCAL MESH (BATHYMETRIC DEPTH)                                      |
  !==============================================================================|
  SUBROUTINE setup_depth
    IMPLICIT NONE
    INTEGER :: IERR
    REAL(SP):: SBUF
    
    IF(dbg_set(dbg_sbr)) write(ipt,*) "SETUP_DEPTH: START"
    hmax = maxval(h_lcl(1:mt))
    hmin = minval(h_lcl(1:mt))


    h = h_lcl

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

  END SUBROUTINE setup_depth

!==============================================================================|
!   SET UP LOCAL CORIOLIS FORCE                                                |
!==============================================================================|
  SUBROUTINE setup_coriolis 
    IMPLICIT NONE
    integer:: I
    !--------------READ IN CORIOLIS PARAMETER--------------------------------------!

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

    CALL n2e2d(c_lcl,cor)


    cor = 2.*7.292e-5_sp * sin(cor * deg2rad)
    
    !!  ggao for equatoral min (4deg)
    IF(.NOT. equator_beta_plane)THEN
     WHERE(cor <  1.e-5_sp .AND. cor > 0.0_sp) cor =  1.e-5_sp
     WHERE(cor > -1.e-5_sp .AND. cor < 0.0_sp) cor = -1.e-5_sp
    END IF 

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

  END SUBROUTINE setup_coriolis
    
    
!==============================================================================|
!   COMPUTE GRAVITY VARIED WITH LATITUDE                                       |
!==============================================================================|

  SUBROUTINE setup_gravity
    IMPLICIT NONE
    INTEGER I

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

    grav_n = grav
    grav_e = grav

    IF(dbg_set(dbg_sbr)) write(ipt,*) "SETUP_GRAVITY: END"
  END SUBROUTINE setup_gravity
    
!==============================================================================|
!   COMPUTE SPONGE LAYER FOR OPEN BOUNDARY DAMPING                             |
!==============================================================================|

  SUBROUTINE setup_sponge
    USE mod_spherical
    IMPLICIT NONE
    REAL(SP)  TEMP,DTMP,C_SPONGE
    INTEGER :: I1, I, SENDER, IERR
    REAL(DP) X1_DP,Y1_DP,X2_DP,Y2_DP,DTMP_DP

    !--SET SPONGE PARAMETERS-------------------------------------------------------|
    
    IF(dbg_set(dbg_sbr)) write(ipt,*) "SETUP_SPONGE: START"

    IF (nsponge ==0) RETURN

! NOTE: X_SPG/Y_SPG COORDINATES MUST BE AJUSTED FOR VXMIN/VYMIN

    x_spg = x_spg - vxmin
    y_spg = y_spg - vymin

    DO i=1,nt
       DO i1=1,nsponge
          dtmp=(xc(i)-x_spg(i1))**2+(yc(i)-y_spg(i1))**2
          dtmp=sqrt(dtmp)/r_spg(i1)

          IF(dtmp <= 1.0_sp) THEN
             c_sponge=c_spg(i1)*(1.0_sp-dtmp)
             cc_sponge(i)=max(c_sponge,cc_sponge(i))
          END IF
       END DO
    END DO


    DEALLOCATE(n_spg,r_spg,c_spg,x_spg,y_spg)

    IF(dbg_set(dbg_sbr)) write(ipt,*) "SETUP_SPONGE: END"
    
    RETURN
  END SUBROUTINE setup_sponge
!==============================================================================|
  
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|

!==============================================================================|
  SUBROUTINE coordinate_units(XL,YL)
    IMPLICIT NONE
    REAL(SP), ALLOCATABLE :: XL(:),YL(:)
    integer status, ierr


    REAL(SP), allocatable :: x_buf(:), lon_buf(:)
    REAL(SP), allocatable :: y_buf(:), lat_buf(:)

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

!! IF NOT SPHERICAL CASE
    
    IF (grid_file_units == 'degrees') THEN
       
       IF (use_proj) THEN
          
          lon = xl
          lat = yl
          
          ! USE PROJECTION TOOL BOX TO CONVERT TO METERS
        IF (serial) THEN
           CALL degrees2meters(xl(1:mt),yl(1:mt), &
               & projection_reference,vx(1:mt),vy(1:mt),mt)
        END IF         
          xm = vx
          ym = vy
          
       ELSE
          
          CALL fatal_error('You must specify a valid projection reference'&
               &,'and compile with PROJ to use files with latitude and longitue in cartesian mode')
       END IF
       

    ELSE IF(grid_file_units == 'meters') THEN

       vx = xl
       vy = yl
       
       xm = xl
       ym = yl
       
       ! USE PROJECTION TOOL BOX TO CONVERT TO DEGREES
       IF (use_proj ) THEN
        IF (serial) THEN
           CALL meters2degrees(xl(1:mt),yl(1:mt), &
            & projection_reference,lon(1:mt),lat(1:mt),mt)
        END IF
          
       END IF

    ELSE 
       CALL fatal_error('UNRECOGNIZED GRID_FILE_UNITS: '//trim(grid_file_units))
       
    END IF


    IF(dbg_set(dbg_sbr)) write(ipt,*) "COODINATE_UNITS: END"
    
  END SUBROUTINE coordinate_units

!==============================================================================|
! SETUP THE SIGMA COORDINATES FOR THE MODEL                                    |
!==============================================================================|
!==============================================================================|
! This program is used to set up the coordinate in the vertical.               !
!                                              !
! case(1): sigma levels                                                        !
! sigma levels are determined by a formula of                                  !
!                      sigma(k)=-[(k-1)/(kb-1)]^k11                            !
!    p_sigma=1: uniform sigma layers                                           !
!    p_sigma=2: layers satisfying a parabolic function with high               !
!               vertical resolution near the surface and bottom.               !
!    p_sigma can be used any real number                                       !
!                                          !
! case(2): general vertical level                                              !
! vertical levels are determined by the formula                                !
!                tanh[(dl+du)((kbm1-k)/kbm1)-dl]+tanh(dl)                      !
!        z(k)= ------------------------------------------  - 1                 !
!                      tanh(dl) + tanh(du)                                     !
!                                                                              !
! case(3): constant layer transformation                                       !
! four values need to be specified:                                            !
!  DUU the upper boundaries, up to which the co-ordinates are parallel must be !
!      defined.                                                                !
!  DLL the lower boundaries, up to which the co-ordinates are parallel must be !
!      defined.                                                                !
!  HMIN1 the minimum water depth at which the layers are constant. If H < HMIN1!
!      then sigma co-ordinates are used.                                       !
!                                                                              !
! Reference of case(2), case(3) and case(4):                                   !
! Pietrzak, J.D., Jan B. Jakobson, Hans Burchard, Hans Jacob Vested, Ole       !
! Petersen , 2002. A three-dimensional hydrostatic model for coastal and ocean !
! modelling using a generalised topography following co-ordinate system. Ocean !
! Modelling 4, 173-205                                                         !
!                                                                              !
!  calculates: z(m,kb) vertical levels                         !
!  calculates: dz(m,kb-1) delta between vertical levels                    !
!  calculates: zz(m,kb-1) intra-vertical levels                    !
!  calculates: dzz(m,kb-2) delta between intra-vertical levels             !
!==============================================================================|

  SUBROUTINE setup_sigma          
    !==============================================================================|
    IMPLICIT NONE
    INTEGER :: K,KK
    INTEGER :: I
    REAL(SP):: ZTMP(KB)
    REAL(SP):: X1,X2,X3  
    REAL(SP):: DR,RCL,RCU
    !==============================================================================|

    IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"SETUP_SIGMA: START"

    IF(dbg_set(dbg_sbrio)) THEN
       WRITE(ipt,*)"==================="
       WRITE(ipt,*)"  SET_SIGMA IO"
       WRITE(ipt,*)"  STYPE = "//trim(stype)
       WRITE(ipt,*)"  P_SIGMA = ", p_sigma
       WRITE(ipt,*)"  KB = ", kb
       WRITE(ipt,*)"==================="
    END IF

    !---------------------------------------------
    !---------------------------------------------
    SELECT CASE(stype)
    !---------------------------------------------
    !SIGMA_COORDINATE_TYPE = UNIFORM (DEGENERATE CASE OF GEOMETRIC)
    CASE(stype_uniform)
    !---------------------------------------------
       IF(p_sigma > 1 .AND. mod(kb,2) == 0) &
            CALL fatal_error('SETUP_SIGMA: COORDINATE TYPE:'//trim(stype)&
            &,'kb shoude be an odd number for this type of sigma coordinates....' )
       CALL sigma_geometric
    !---------------------------------------------
    !SIGMA_COORDINATE_TYPE = GEOMETRIC
    CASE(stype_geometric)
    !---------------------------------------------

       IF(p_sigma > 1 .AND. mod(kb,2) == 0) &
            CALL fatal_error('SETUP_SIGMA: COORDINATE TYPE:'//trim(stype)&
            &,'kb shoude be an odd number for this type of sigma coordinates....' )
       CALL sigma_geometric
    !---------------------------------------------
    ! SIGMA_COORDINATE_TYPE = TANH
    CASE(stype_tanh)
    !---------------------------------------------   
       CALL sigma_tanh
    !---------------------------------------------   
    !SIGMA_COORDINATE_TYPE = GENERALIZED
    CASE(stype_generalized)
    !---------------------------------------------   
       CALL sigma_generalized

    ! THIS IS A CURRENTLY UNUSED METHOD
!!$!---------------------------------------------   
!!$!SIGMA_COORDINATE_TYPE = WHAT THE HELL IS THIS?
!!$   CASE("UNKNOWN")
!!$!---------------------------------------------   
!!$
!!$    CALL SIGMA_UNKNOWN
    CASE DEFAULT
       CALL fatal_error("SET_SIGMA: REACHED DEFAULT CASE FOR SIGMA COOR&
            &DINATE TYPE")
    END SELECT

    !---------COMPUTE SIGMA DERIVATIVES AND INTRA SIGMA LEVELS---------------------!

    IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"SETUP_SIGMA: END"

    RETURN

  END SUBROUTINE setup_sigma
!==============================================================================|
  SUBROUTINE sigma_geometric
      IMPLICIT NONE
      INTEGER :: I,K
      REAL(SP):: ZTMP(KB)
      !orginal formula to set sigma   
      IF(p_sigma == 1)THEN
         DO k=1,kb
            ztmp(k) = -((k-1)/float(kb-1))**p_sigma 
         END DO
      ELSE
         DO k=1,(kb+1)/2
            ztmp(k) = -((k-1)/float((kb+1)/2-1))**p_sigma/2 
         END DO
         DO k=(kb+1)/2+1,kb
            ztmp(k) = ((kb-k)/float((kb+1)/2-1))**p_sigma/2-1.0
         END DO
      END IF

      DO i=1,m
         DO k=1,kb
            z(i,k)=ztmp(k)
         END DO
      END DO

      DO i=1,n
         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
    END SUBROUTINE sigma_geometric
!--------------------------------------------------------------------    
    SUBROUTINE sigma_generalized
      IMPLICIT NONE
      INTEGER :: I,K, kk
      REAL(SP):: X1,X2,X3  
      REAL(SP):: DR,RCL,RCU
 
      DO i=1,m
         IF(h(i) < hmin1)THEN
            z(i,1)=0.0
            dl2=0.001;du2=0.001
            DO k=1,kbm1
               x1=dl2+du2
               x1=x1*(kbm1-k)/kbm1
               x1=x1-dl2
               x1=tanh(x1)
               x2=tanh(dl2)
               x3=x2+tanh(du2)

               z(i,k+1)=(x1+x2)/x3-1.0_sp
            END DO
         ELSE
            dr=(h(i)-duu-dll)/h(i)/(kb-ku-kl-1)

            z(i,1)=0.0_sp

            DO k=2,ku+1
               z(i,k)=z(i,k-1)-zku(k-1)/h(i)
            END DO

            DO k=ku+2,kb-kl
               z(i,k)=z(i,k-1)-dr
            END DO

            kk=0
            DO k=kb-kl+1,kb
               kk=kk+1
               z(i,k)=z(i,k-1)-zkl(kk)/h(i)
            END DO
         END IF
      END DO

      DO i=1,n
         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
    END SUBROUTINE sigma_generalized
!--------------------------------------------------------------------
    SUBROUTINE sigma_tanh
      IMPLICIT NONE
      INTEGER :: I,K
      REAL(SP):: X1,X2,X3  
    
      z=0.0;z1=0.0
      DO k=1,kbm1
         x1=dl2+du2
         x1=x1*(kbm1-k)/kbm1
         x1=x1-dl2
         x1=tanh(x1)
         x2=tanh(dl2)
         x3=x2+tanh(du2)
         DO i=1,m
            z(i,k+1)=(x1+x2)/x3-1.0_sp
         END DO
         DO i=1,n
            z1(i,k+1)=(x1+x2)/x3-1.0_sp
         END DO
      END DO
    END SUBROUTINE sigma_tanh
!--------------------------------------------------------------------
!!$ SUBROUTINE SIGMA_UNKNOWN
!!$   DO I=1,M
!!$     IF(H(I) < HMIN1)THEN
!!$       RCU=-DUU/HMIN1
!!$       RCL=DLL/HMIN1-1
!!$       DR=(RCL-RCU)/(KB-KU-KL-1)
!!$
!!$       DO K=1,KU
!!$         ZKU(K)=RCU/KU
!!$       END DO
!!$       DO K=1,KL
!!$         ZKL(K)=(-1.0_SP-RCL)/KL
!!$       END DO
!!$
!!$       Z(I,1)=0.0_SP
!!$
!!$       DO K=2,KU+1
!!$         Z(I,K)=Z(I,K-1)+ZKU(K-1)
!!$       END DO
!!$
!!$       DO K=KU+2,KB-KL
!!$         Z(I,K)=Z(I,K-1)+DR
!!$       END DO
!!$
!!$       KK=0
!!$       DO K=KB-KL+1,KB
!!$         KK=KK+1
!!$         Z(I,K)=Z(I,K-1)+ZKL(KK)
!!$       END DO
!!$
!!$     ELSE
!!$       DR=(H(I)-DUU-DLL)/H(I)/(KB-KU-KL-1)
!!$
!!$       Z(I,1)=0.0_SP
!!$
!!$       DO K=2,KU+1
!!$         Z(I,K)=Z(I,K-1)-ZKU(K-1)/H(I)
!!$       END DO
!!$
!!$       DO K=KU+2,KB-KL
!!$         Z(I,K)=Z(I,K-1)-DR
!!$       END DO
!!$
!!$       KK=0
!!$       DO K=KB-KL+1,KB
!!$         KK=KK+1
!!$         Z(I,K)=Z(I,K-1)-ZKL(KK)/H(I)
!!$       END DO
!!$     END IF
!!$   END DO
!!$   
!!$   DO I=1,N
!!$     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    
!!$END SUBROUTINE SIGMA_UNKNOWN

!==============================================================================|

  SUBROUTINE setup_sigma_derivatives
    USE all_vars
    IMPLICIT NONE
    INTEGER :: K, I


    IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"SETUP_SIGMA_DERIVATIVES: START"

    DO k=1,kb-1
       DO i=1,mt
          dz(i,k)  = z(i,k)-z(i,k+1)
          zz(i,k)  = .5_sp*(z(i,k)+z(i,k+1))
       END DO
       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,mt
       zz(i,kb) = 2.0_sp*zz(i,kb-1)-zz(i,kb-2)
    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,mt
          dzz(i,k) = zz(i,k)-zz(i,k+1)
       END DO
       DO i=1,nt
          dzz1(i,k) = zz1(i,k)-zz1(i,k+1)
       END DO
    END DO

    dzz(:,kbm1) = 0.0_sp
    dz(:,kb)    = 0.0_sp
    dzz1(:,kbm1) = 0.0_sp
    dz1(:,kb)    = 0.0_sp


    !----------OUTPUT VALUES-TO INFOFILE-------------------------------------------!

    IF(dbg_set(dbg_log)) THEN
       WRITE(ipt,*  )'!'
       WRITE(ipt,*  )'!'
       WRITE(ipt,*)'!                SIGMA LAYER INFO     '
       WRITE(ipt,*) "SIGMA TYPE:",trim(stype)
       WRITE(ipt,70)
       SELECT CASE(stype)
       CASE(stype_uniform)
          DO k=1,kb
             WRITE(ipt,80) k,z(1,k),zz(1,k),dz(1,k),dzz(1,k)
          END DO
       CASE(stype_restart)
          DO k=1,kb
             WRITE(ipt,80) k,z(1,k),zz(1,k),dz(1,k),dzz(1,k)
          END DO
       CASE(stype_geometric)
          DO k=1,kb
             WRITE(ipt,80) k,z(1,k),zz(1,k),dz(1,k),dzz(1,k)
          END DO
       CASE(stype_tanh)
          DO k=1,kb
             WRITE(ipt,80) k,z(1,k),zz(1,k),dz(1,k),dzz(1,k)
          END DO
       CASE(stype_generalized) ! THIS IS CASE SPECIFIC 
          WRITE(ipt,*)"SET CASE SPECIFIC GENERALIZED SIGMA LAYER OUTPUT TO SCREEN &
               &IN mod_setup.F"
       END SELECT
       WRITE(ipt,*  )'!'
    END IF

    IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"END SETUP_SIGMA_DERIVATIVES"

    !----------FORMAT STATEMENTS---------------------------------------------------!

70  FORMAT(2x,'k',13x,'z',11x,'zz',11x,'dz',11x,'dzz')
80  FORMAT(' ',i5,4f13.8)


    IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"SETUP_SIGMA_DERIVATIVES: END"
  END SUBROUTINE setup_sigma_derivatives
  
END MODULE mod_setup