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

!==============================================================================!
! TYPE_OBC = 1: Surface Elevation Specified (Tidal Forcing) (ASL)              !
! TYPE_OBC = 2: AS TYPE_OBC=1 AND NON-LINEAR FLUX FOR CURRENT AT OPEN BOUNDARY !
! TYPE_OBC = 3: Zero Surface Elevation Boundary Condition (ASL_CLP)            !
! TYPE_OBC = 4: AS TYPE_OBC=3 AND NON-LINEAR FLUX FOR CURRENT AT OPEN BOUNDARY !
! TYPE_OBC = 5: GRAVITY-WAVE RADIATION IMPLICIT OPEN BOUNDARY CONDITION (GWI)  !
! TYPE_OBC = 6: AS TYPE_OBC=5 AND NON-LINEAR FLUX FOR CURRENT AT OPEN BOUNDARY !
! TYPE_OBC = 7: BLUMBERG AND KHANTA IMPLICIT OPEN BOUNDARY CONDITION (BKI)     !
! TYPE_OBC = 8: AS TYPE_OBC=7 AND NON-LINEAR FLUX FOR CURRENT AT OPEN BOUNDARY !
! TYPE_OBC = 9: ORLANSKI RADIATION EXPLICIT OPEN BOUNDARY CONDITION (ORE)      !
! TYPE_OBC =10: AS TYPE_OBC=9 AND NON-LINEAR FLUX FOR CURRENT AT OPEN BOUNDARY !
!                                                                              !
! TYPE_TSOBC = 1: THE PERTURBATION OF TEMPERATURE AND SALINITY AT OPEN BOUNDARY!
!                 ARE EQUAL TO THAT AT NEXT_OBC                                !
! TYPE_TSOBC = 2: GRAVITY-WAVE RADIATION IMPLICIT OPEN BOUNDARY CONDITION FOR  !
!                 THE PERTURBATION OF TEMPERATURE AND SALINITY                 !
! TYPE_TSOBC = 3: BLUMBERG AND KHANTA IMPLICIT OPEN BOUNDARY CONDITION FOR     !
!                 THE PERTURBATION OF TEMPERATURE AND SALINITY                 !
! TYPE_TSOBC = 4: ORLANSKI RADIATION EXPLICIT OPEN BOUNDARY CONDITION FOR      !
!                 THE PERTURBATION OF TEMPERATURE AND SALINITY                 !
!==============================================================================!

MODULE mod_obcs
  USE mod_prec
  USE mod_utils
  USE control,  type_tsobc => obc_ts_type
  USE mod_time
  USE mod_force
  USE bcs
  IMPLICIT NONE
  SAVE

  !--Open Boundary Types, Lists, Pointers

  !   INTEGER               :: IOBCN_GL         !!GLOBAL NUMBER OF OPEN BOUNDARY NODES
  !   INTEGER               :: IOBCN            !!LOCAL NUMBER OF OPEN BOUNDARY NODES
  !   INTEGER,  ALLOCATABLE :: I_OBC_GL(:)      !!GLOBAL ID OF OPEN BOUNDARY NODES
  !   INTEGER,  ALLOCATABLE :: I_OBC_N(:)       !!OPEN BOUNDARY NODE LIST
  INTEGER,  ALLOCATABLE :: next_obc(:)      !!INTERIOR NEIGHBOR OF OPEN BOUNDARY NODE
  INTEGER,  ALLOCATABLE :: next_obc2(:)     !!INTERIOR NEIGHBOR OF NEXT_OBC
  !   INTEGER,  ALLOCATABLE :: TYPE_OBC(:)      !!OUTER BOUNDARY NODE TYPE (FOR SURFACE ELEVATION)
  !   INTEGER,  ALLOCATABLE :: TYPE_OBC_GL(:)   !!OUTER BOUNDARY NODE TYPE (FOR SURFACE ELEVATION)
  INTEGER               :: ibcn(5)          !!NUMBER OF EACH TYPE OF OBN IN LOCAL  DOM
  INTEGER               :: ibcn_gl(5)       !!NUMBER OF EACH TYPE OF OBN IN GLOBAL DOM
  INTEGER,  ALLOCATABLE :: obc_lst(:,:)     !!MAPPING OF OPEN BOUNDARY ARRAYS TO EACH TYPE
  INTEGER,  ALLOCATABLE :: nadjn_obc(:)     !!NUMBER OF ADJACENT OPEN BOUNDARY NODES TO OBN
  INTEGER,  ALLOCATABLE :: adjn_obc(:,:)    !!ADJACENT OBN's of OBN
  INTEGER,  ALLOCATABLE :: nadjc_obc(:)     !!NUMBER OF ADJACENT OPEN BOUNDARY CELLS TO OBN
  INTEGER,  ALLOCATABLE :: adjc_obc(:,:)    !!ADJACENT OPEN BOUNDARY CELLS

  !-- Open Boundary TEMPERATURE AND SALINITY FOR NUDGING
  REAL(sp), ALLOCATABLE :: temp_obc(:,:)      !!OPEN BOUNDARY TEMPERATURE (see mod_force.F)     
  REAL(sp), ALLOCATABLE :: salt_obc(:,:)      !!OPEN BOUNDARY SALT (see mod_force.F)   
  !--Open Boundary Metrics
  INTEGER,  ALLOCATABLE :: nfluxf_obc(:)    !!NUMBER OF FLUX SEGMENTS TO OBN
  REAL(sp), ALLOCATABLE :: fluxf_obc(:,:)   !!FLUX FRACTION ON EACH SIDE OF OBN
  REAL(sp), ALLOCATABLE :: nxobc(:)         !!INWARD POINTING X-NORMAL OF OBN
  REAL(sp), ALLOCATABLE :: nyobc(:)         !!INWARD POINTING Y-NORMAL OF OBN
  REAL(sp), ALLOCATABLE :: dltn_obc(:)      !!DISTANCE BETWEEN NEXT_OBC AND OBN NORMAL TO BOUNDARY

  !--Previous Time Level Free Surface Fields
  REAL(sp), ALLOCATABLE :: elm1(:)          !!SURFACE ELEV FROM PREVIOUS TIME LEVEL (ORLANSKI COND)
  REAL(sp), ALLOCATABLE :: elm2(:)          !!SURFACE ELEV FROM PREVIOUS TIME LEVEL (ORLANSKI COND)
  REAL(sp), ALLOCATABLE :: t1m1(:,:)          !!TEMPERATURE FROM PREVIOUS TIME LEVEL (ORLANSKI COND)
  REAL(sp), ALLOCATABLE :: t1m2(:,:)          !!TEMPERATURE FROM PREVIOUS TIME LEVEL (ORLANSKI COND)
  REAL(sp), ALLOCATABLE :: s1m1(:,:)          !!SALINITY FROM PREVIOUS TIME LEVEL (ORLANSKI COND)
  REAL(sp), ALLOCATABLE :: s1m2(:,:)          !!SALINITY FROM PREVIOUS TIME LEVEL (ORLANSKI COND)

  !--Nonlinear Velocity Open Boundary Condition Arrays
  REAL(sp), ALLOCATABLE :: fluxobn(:)
  REAL(sp), ALLOCATABLE :: iucp(:)
  REAL(sp), ALLOCATABLE :: xflux_obcn(:)
  REAL(sp), ALLOCATABLE :: uard_obcn(:)
  REAL(sp), ALLOCATABLE :: xflux_obc(:,:)


CONTAINS
  !==========================================================================|
  SUBROUTINE alloc_obc_data

    !--------------------------------------------------------------------------|
    !  ALLOCATE AND INITIALIZE SURFACE ELEVATION ARRAYS FOR                    |
    !  TIME STEPS N-1 AND N-2                                                  |
    !--------------------------------------------------------------------------|

    USE all_vars
    IMPLICIT NONE

    ALLOCATE(elm1(0:mt))             ;elm1       = zero 
    ALLOCATE(elm2(0:mt))             ;elm2       = zero 
    ALLOCATE(next_obc(iobcn))        ;next_obc   = 0
    ALLOCATE(next_obc2(iobcn))       ;next_obc2  = 0
    ALLOCATE(fluxobn(1:nt))          ;fluxobn    = zero
    ALLOCATE(iucp(0:nt))             ;iucp       = 1
    ALLOCATE(xflux_obcn(iobcn+1))    ;xflux_obcn = zero
    ALLOCATE(uard_obcn(iobcn+1))     ;uard_obcn  = zero
    ALLOCATE(xflux_obc(iobcn+1,kbm1));xflux_obc  = zero

    ALLOCATE(temp_obc(iobcn,kbm1))   ;temp_obc   = zero
    ALLOCATE(salt_obc(iobcn,kbm1))   ;salt_obc   = zero

    IF(type_tsobc == 4)THEN
       ALLOCATE(t1m1(0:mt,kbm1))        ;t1m1       = zero 
       ALLOCATE(t1m2(0:mt,kbm1))        ;t1m2       = zero 
       ALLOCATE(s1m1(0:mt,kbm1))        ;s1m1       = zero 
       ALLOCATE(s1m2(0:mt,kbm1))        ;s1m2       = zero 
    END IF
    RETURN
  END SUBROUTINE alloc_obc_data
  !==========================================================================|

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

  !==========================================================================|
  SUBROUTINE assign_elm1_to_elm2

    !--------------------------------------------------------------------------|
    !  Assign ELM1 to ELM2 and EL to ELM1                                      |
    !--------------------------------------------------------------------------|

    USE all_vars
    IMPLICIT NONE

    elm2 = elm1
    elm1 = el

    IF(type_tsobc == 4)THEN
       t1m2 = t1m1
       t1m1 = t1

       s1m2 = s1m1
       s1m1 = s1
    END IF

    RETURN
  END SUBROUTINE assign_elm1_to_elm2
  !==========================================================================|

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


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

  !==========================================================================|
  SUBROUTINE elevation_atmo ! ONLY USED if defined (ATMO_TIDE)
    !--------------------------------------------------------------------------|
    !  Surface Elevation of ATMOSPHERIC TIDE                                   |
    !--------------------------------------------------------------------------|
    USE all_vars
    IMPLICIT NONE

    REAL(SP),PARAMETER :: APT_ATMO = 0.0113_sp       
    REAL(DP),PARAMETER :: ALFA_ATMO = 112.0_sp*pi2/360.0_sp
    REAL(SP),PARAMETER :: S2PERIOD = 43200.0_sp

    INTEGER :: I,J
    TYPE(time):: TIME_ELAPSED
    REAL(DP):: TIME_FLT,PHAI_IJ
    REAL(SP):: FORCE

    IF(.not.use_proj) CALL fatal_error &
         & ("ELEVATION_ATMO: ILLEGAL ATTEMPT TO USE ATMOSPHERIC TIDE IN CARTESIAN MODE",&
         & "YOU MUST BE USING PROJ4 TO DO THIS",&
         & "THE POLICE HAVE BEEN NOTIFIED AND WILL BE THERE SHORTLY")




    ! THIS CODE IS CALLED DURING EACH RK STAGE
    ! DECIDE WHICH TIME YOU WANT TO USE!
    time_elapsed = rktime - spectime

!! IN MY JUDGEMENT, A TIMESERIES FORCING DOES NOT MAKE SENSE FOR ATMO/EQUITIDE

!    SELECT CASE(TIDE_FORCING_TYPE)
!    CASE(TIDE_FORCING_TIMESERIES)
       !
       !-Julian: Set Elevation Based on Linear Interpolation Between Two Data Times-|
       !
!       CALL FATAL_ERROR("ELEVATION_EQUI: Not set up for Julian Timeseries forcing yet?")

!    CASE(TIDE_FORCING_SPECTRAL)
       !
       !-Non-Julian: Set Elevation of Equilibrium Tide -----------------------------|
       !

       time_flt = seconds(time_elapsed * pi2/s2period)

       DO i = 1, mt
          phai_ij = lon(i)*pi2/360.0_sp
          force = apt_atmo*dcos(time_flt+2.0_dp*phai_ij-alfa_atmo)
          elf_atmo(i) = force * ramp
       END DO

!    CASE DEFAULT
!       CALL FATAL_ERROR("ELEVATION_ATMO: INVALID TIDE FORCING TYPE ?")
!    END SELECT



    RETURN
  END SUBROUTINE elevation_atmo
  !==========================================================================|

!==============================================================================|
!     SET BOUNDARY CONDITIONS FOR CALCULATING atmospher pressure external model|
!                                                                              |
!==============================================================================|

   SUBROUTINE bcond_pa_air
   END SUBROUTINE bcond_pa_air
   
   
!==============================================================================|
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|

  !==========================================================================|
  SUBROUTINE bcond_asl

    !--------------------------------------------------------------------------|
    !  Surface Elevation Boundary Condition (Tidal Forcing)                    |
    !--------------------------------------------------------------------------|

    USE all_vars, only : elf
    IMPLICIT NONE

    INTEGER :: I,II,J,JN
    TYPE(time):: TIME_ELAPSED
    REAL(DP):: TIME_FLT,PHAI_IJ
    REAL(SP):: FORCE
    REAL(SP), ALLOCATABLE :: BND_ELV(:)


    ALLOCATE(bnd_elv(iobcn))
    !
    !-Julian: Set Elevation Based on Linear Interpolation Between Two Data Times-|
    !
    SELECT CASE(tide_forcing_type)
    CASE(tide_forcing_timeseries)


       ! THIS CODE IS CALLED DURING EACH RK STAGE
       ! DECIDE WHICH TIME YOU WANT TO USE!
       CALL update_tide(rktime,bnd_elv)
       !      CALL UPDATE_TIDE(ExtTime,BND_ELV)


       DO j = 1, ibcn(1)
          jn = obc_lst(1,j)
          ii = i_obc_n(jn)

          elf(ii) = bnd_elv(j)*ramp
       END DO

       !
       !-Non-Julian: Set Elevation Based on Input Amplitude and Phase of Tidal Comps-|
       !

    CASE(tide_forcing_spectral)




       ! THIS CODE IS CALLED DURING EACH RK STAGE
       ! DECIDE WHICH TIME YOU WANT TO USE!
       time_elapsed = rktime - spectime
       !    TIME_ELAPSED = ExtTime - SpecTime

       !      write(ipt,*) "-----------------------------------------"
       DO i = 1, ibcn(1)
          jn = obc_lst(1,i)
          ii = i_obc_n(jn)
          force = 0.0_sp
          DO j = 1,ntidecomps
             phai_ij = phai(i,j)*pi2/360.0_sp
             time_flt = seconds(time_elapsed * pi2/period(j))
             ! Take cosine of a double precision number to preserve accuracy!
             force = apt(i,j)*dcos(time_flt -phai_ij) + force
          END DO
          force = force + emean(i)
          elf(ii) = force * ramp
          !         write(ipt,*) "ELF=",ELF(II)

       END DO

    CASE DEFAULT
       CALL fatal_error("BCOND_ASL: INVALID TIDE FORCING TYPE ?")
    END SELECT

    DEALLOCATE(bnd_elv)

    RETURN
  END SUBROUTINE bcond_asl
!!$!==========================================================================|

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

  !==========================================================================|
  SUBROUTINE bcond_asl_clp

    !--------------------------------------------------------------------------|
    !  Zero Surface Elevation Boundary Condition                               |
    !--------------------------------------------------------------------------|

    USE all_vars
    IMPLICIT NONE

    INTEGER :: I,II,JN

    DO i = 1, ibcn(2)
       jn = obc_lst(2,i)
       ii = i_obc_n(jn)  
       elf(ii) = 0.0_sp
    END DO

    RETURN
  END SUBROUTINE bcond_asl_clp
  !==========================================================================|

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

  !==========================================================================|
  SUBROUTINE bcond_gwi(K_RK)

    !--------------------------------------------------------------------------|
    !  GRAVITY-WAVE RADIATION IMPLICIT OPEN BOUNDARY CONDITION (GWI)           |
    !--------------------------------------------------------------------------|

    USE all_vars
    IMPLICIT NONE

    INTEGER :: I1,I2,J,JN,K_RK
    REAL(SP):: CC,CP,ALPHA_RK_TMP

    IF(k_rk == 0)THEN
      alpha_rk_tmp = 1.0
    ELSE
      alpha_rk_tmp = alpha_rk(k_rk)
    END IF

    DO j = 1,ibcn(3)
      jn = obc_lst(3,j)
      i1 =  i_obc_n(jn) 
      i2 = next_obc(jn) 
      cc = sqrt(grav_n(i1)*h(i1))*dte/dltn_obc(jn)*alpha_rk_tmp
      cp = cc + 1.0_sp
      elf(i1) = (cc*elf(i2) + elrk(i1))/cp
    END DO

    RETURN
  END SUBROUTINE bcond_gwi
  !==========================================================================|

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

  !==========================================================================|
  SUBROUTINE bcond_bki(K_RK)

    !--------------------------------------------------------------------------|
    !  BLUMBERG AND KHANTA IMPLICIT OPEN BOUNDARY CONDITION                    |
    !--------------------------------------------------------------------------|

    USE all_vars
    IMPLICIT NONE

    INTEGER :: I1,I2,J,JN,K_RK
    REAL(SP):: CC,CP,ALPHA_RK_TMP

    IF(k_rk == 0)THEN
      alpha_rk_tmp = 1.0
    ELSE
      alpha_rk_tmp = alpha_rk(k_rk)
    END IF

    DO j = 1,ibcn(4)
      jn = obc_lst(4,j)
      i1 = i_obc_n(jn)
      i2 = next_obc(jn)
      cc = sqrt(grav_n(i1)*h(i1))*dte/dltn_obc(jn)*alpha_rk_tmp
      cp = cc + 1.0_sp
      elf(i1) = (cc*elf(i2) + elrk(i1)*(1.0_sp-dte/10800.0_sp))/cp
    END DO

    RETURN
  END SUBROUTINE bcond_bki
  !==============================================================================|

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

  !==============================================================================|
  SUBROUTINE bcond_ore

    !------------------------------------------------------------------------------|
    !  ORLANSKI RADIATION EXPLICIT OPEN BOUNDARY CONDITION (ORE)                   |
    !------------------------------------------------------------------------------|

    USE all_vars
    IMPLICIT NONE

    INTEGER  :: I1,I2,I3,J,JN
    REAL(SP) :: CL, MU

    DO j = 1, ibcn(5)
       jn = obc_lst(5,j)
       i1 = i_obc_n(jn)
       i2 = next_obc(jn)
       i3 = next_obc2(jn)

       cl = (elm2(i2)-el(i2))/(el(i2)+elm2(i2)-2.0*elm1(i3))
       IF(cl >= 1.0)THEN
          mu = 1.0
       ELSE IF(cl > 0.0 .AND. cl < 1.0)THEN
          mu = cl
       ELSE
          mu = 0.0
       END IF

       elf(i1)=(elm1(i1)*(1.0-mu)+2.0*mu*el(i2))/(1.0+mu)
    END DO

    RETURN
  END SUBROUTINE bcond_ore
  !==============================================================================|

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

  !========================================================================

  SUBROUTINE setup_obctypes
    USE lims
    IMPLICIT NONE

    INTEGER              :: I,I1,I2,NCNT,IERR,J
    INTEGER, ALLOCATABLE :: TEMP1(:),TEMP2(:),TEMP3(:),TEMP4(:),TEMP5(:),TEMP6(:),TEMP7(:),ITEMP(:)

    ! CREATE THE BC TYPE INTEGERS USED IN THE MODEL
    CALL separate_obc ! (MOD_OBCS.F)

    !==============================================================================|
    !   REPORT AND CHECK OBC RESULTS                                               |
    !==============================================================================|
    ALLOCATE(temp1(nprocs),temp2(nprocs),temp3(nprocs),temp4(nprocs))
    ALLOCATE(temp5(nprocs),temp6(nprocs),temp7(nprocs))
    temp1(1)  = iobcn
    temp2(1) = ibcn(1)
    temp3(1) = ibcn(2)
    temp4(1) = ibcn(3)
    temp5(1) = ibcn(4)
    temp6(1) = ibcn(5)
    !   TEMP7(1) = NOBCGEO



    IF(dbg_set(dbg_log)) THEN
       WRITE(ipt,*)  '! BCOND TYPE             : TOTAL BC NODES; BC NODES IN EACH DOMAIN'
       WRITE(ipt,100)'!  IOBCN                 :',iobcn_gl,   (temp1(i),i=1,nprocs)
       WRITE(ipt,100)'!  IBCN(1)               :',ibcn_gl(1), (temp2(i),i=1,nprocs)
       WRITE(ipt,100)'!  IBCN(2)               :',ibcn_gl(2), (temp3(i),i=1,nprocs)
       WRITE(ipt,100)'!  IBCN(3)               :',ibcn_gl(3), (temp4(i),i=1,nprocs)
       WRITE(ipt,100)'!  IBCN(4)               :',ibcn_gl(4), (temp5(i),i=1,nprocs)
       WRITE(ipt,100)'!  IBCN(5)               :',ibcn_gl(5), (temp6(i),i=1,nprocs)
       !   IF(MSR)WRITE(IPT,100)'!  NOBCGEO               :'
       !,NOBCGEO_GL, (TEMP7(I),I=1,NPROCS)
    END IF

    DEALLOCATE(temp1,temp2,temp3,temp4,temp5,temp6,temp7)

    RETURN
100 FORMAT(1x,a26,i6," =>",2x,4(i5,1h,))

  END SUBROUTINE setup_obctypes
  !==============================================================================|

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

  !==============================================================================|
  SUBROUTINE separate_obc

    !------------------------------------------------------------------------------|
    ! Accumulate separately the amounts of nodes for 11 types of open boundaries   |
    !------------------------------------------------------------------------------|
    USE all_vars
    IMPLICIT NONE

    INTEGER :: I,I1,I2,I3,I4,I5,II,J

    ibcn = 0
    ibcn_gl = 0

    DO i = 1, iobcn_gl
       IF(type_obc_gl(i) == 1 .OR. type_obc_gl(i) == 2)  ibcn_gl(1) = ibcn_gl(1) + 1
       IF(type_obc_gl(i) == 3 .OR. type_obc_gl(i) == 4)  ibcn_gl(2) = ibcn_gl(2) + 1
       IF(type_obc_gl(i) == 5 .OR. type_obc_gl(i) == 6)  ibcn_gl(3) = ibcn_gl(3) + 1
       IF(type_obc_gl(i) == 7 .OR. type_obc_gl(i) == 8)  ibcn_gl(4) = ibcn_gl(4) + 1
       IF(type_obc_gl(i) == 9 .OR. type_obc_gl(i) == 10) ibcn_gl(5) = ibcn_gl(5) + 1
    END DO

    DO i = 1, iobcn
       IF(type_obc(i) == 1 .OR. type_obc(i) == 2)  ibcn(1) = ibcn(1) + 1
       IF(type_obc(i) == 3 .OR. type_obc(i) == 4)  ibcn(2) = ibcn(2) + 1
       IF(type_obc(i) == 5 .OR. type_obc(i) == 6)  ibcn(3) = ibcn(3) + 1
       IF(type_obc(i) == 7 .OR. type_obc(i) == 8)  ibcn(4) = ibcn(4) + 1
       IF(type_obc(i) == 9 .OR. type_obc(i) == 10) ibcn(5) = ibcn(5) + 1
    END DO

    i1 = 0
    i2 = 0
    i3 = 0
    i4 = 0
    i5 = 0


    ALLOCATE(obc_lst(5,maxval(ibcn))) ; obc_lst = 0

    DO i=1,iobcn
       IF(type_obc(i) == 1 .OR. type_obc(i) == 2)THEN
          i1 = i1 + 1
          obc_lst(1,i1) = i
       ELSE IF(type_obc(i) == 3 .OR. type_obc(i) == 4)THEN
          i2 = i2 + 1
          obc_lst(2,i2) = i
       ELSE IF(type_obc(i) == 5 .OR. type_obc(i) == 6)THEN
          i3 = i3 + 1
          obc_lst(3,i3) = i
       ELSE IF(type_obc(i) == 7 .OR. type_obc(i) == 8)THEN
          i4 = i4 + 1
          obc_lst(4,i4) = i
       ELSE IF(type_obc(i) == 9 .OR. type_obc(i) == 10)THEN
          i5 = i5 + 1
          obc_lst(5,i5) = i
       END IF
    END DO

    RETURN
  END SUBROUTINE separate_obc
  !==============================================================================|

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

  !==============================================================================|
  SUBROUTINE setup_obc

    !------------------------------------------------------------------------------!
    USE all_vars
    IMPLICIT NONE

    REAL(SP) :: DXC,DYC,DXN,DYN,CROSS,E1,E2,DOTMAX,DOT,DX,DY,DXN_TMP,DYN_TMP
    INTEGER  :: I,J,JJ,INODE,JNODE,I1,I2,IC,N1,N2,N3
    LOGICAL  :: DEBUG

    REAL(SP), ALLOCATABLE :: NXOBC_TMP(:),NYOBC_TMP(:)

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


    IF(dbg_set(dbg_log))WRITE(ipt,*) "!"
    IF(dbg_set(dbg_log))WRITE(ipt,*) "!           SETTING UP OPEN BOUNDARY CONDITIONS"
    IF(dbg_set(dbg_log))WRITE(ipt,*) "!"

    !--Determine Adjacent Open Boundary Points-------------------------------------!

    ALLOCATE(nadjn_obc(iobcn))  ; nadjn_obc = 0
    ALLOCATE(adjn_obc(iobcn,2)) ; adjn_obc = 0

    DO i=1,iobcn
       inode = i_obc_n(i)
       DO j=1,ntsn(inode)-1
          jnode = nbsn(inode,j)
          IF(isonb(jnode) == 2 .AND. inode /= jnode)THEN
             nadjn_obc(i) = nadjn_obc(i) + 1
             adjn_obc(i,nadjn_obc(i)) = jnode
          END IF
       END DO
    END DO

    DO i=1,iobcn
       IF(nadjn_obc(i) == 0)THEN
          WRITE(*,*)'NO ADJACENT NODE FOUND FOR BOUNDARY NODE',i
          WRITE(*,*)'IN PROCESSOR',myid
          CALL pstop
       END IF
    END DO

    !--Determine Adjacent Cells-(Nonlinear Only)-----------------------------------!
    !--Simultaneously Determine INWARD Pointing Normal NXOBC,NYOBC                 !

    ALLOCATE(nadjc_obc(iobcn))  ; nadjc_obc = 0
    ALLOCATE(adjc_obc(iobcn,2)) ; adjc_obc = 0
    ALLOCATE(nxobc(iobcn)) ; nxobc = 0
    ALLOCATE(nyobc(iobcn)) ; nyobc = 0
    ALLOCATE(nxobc_tmp(iobcn)) ; nxobc_tmp = 0
    ALLOCATE(nyobc_tmp(iobcn)) ; nyobc_tmp = 0

    DO i=1,iobcn
       i1 = i_obc_n(i)

       !!Mark First Cell on Boundary Edge Adjacent to Node I
       i2 = adjn_obc(i,1)
       DO j = 1, ntve(i1)
          ic = nbve(i1,j)
          n1 = nv(ic,1) ; n2 = nv(ic,2) ; n3 = nv(ic,3)
          IF( n1-i2 == 0 .OR. n2-i2 == 0 .OR. n3-i2 == 0)THEN
             dxn = vx(i2)-vx(i1) ; dyn = vy(i2)-vy(i1)
             dxc = xc(ic)-vx(i1) ; dyc = yc(ic)-vy(i1)
             cross = sign(1.0_sp,dxc*dyn - dyc*dxn)
             nxobc_tmp(i) =  cross*dyn/sqrt(dxn**2 +dyn**2)
             nyobc_tmp(i) = -cross*dxn/sqrt(dxn**2 +dyn**2)
             nxobc(i) = nxobc_tmp(i)
             nyobc(i) = nyobc_tmp(i)
             nadjc_obc(i) = nadjc_obc(i) + 1
             adjc_obc(i,nadjc_obc(i)) = ic
             IF(mod(type_obc(i),2) == 1)THEN
                !!Node is Linear, Mark Cell as Linear for Flux Update
                iucp(ic) = 0
             END IF
          END IF
       END DO

       IF(nadjn_obc(i) > 1)THEN
          i2 = adjn_obc(i,2)
          DO j = 1, ntve(i1)
             ic = nbve(i1,j)
             n1 = nv(ic,1) ; n2 = nv(ic,2) ; n3 = nv(ic,3)
             IF( n1-i2 == 0 .OR. n2-i2 == 0 .OR. n3-i2 == 0)THEN
                dxn = vx(i2)-vx(i1) ; dyn = vy(i2)-vy(i1)
                dxc = xc(ic)-vx(i1) ; dyc = yc(ic)-vy(i1)
                cross = sign(1.0_sp,dxc*dyn - dyc*dxn)
                nxobc_tmp(i) = nxobc_tmp(i) + cross*dyn/sqrt(dxn**2 +dyn**2)
                nyobc_tmp(i) = nyobc_tmp(i) - cross*dxn/sqrt(dxn**2 +dyn**2)
                nxobc(i) = nxobc_tmp(i)/sqrt(nxobc_tmp(i)**2 + nyobc_tmp(i)**2)
                nyobc(i) = nyobc_tmp(i)/sqrt(nxobc_tmp(i)**2 + nyobc_tmp(i)**2)

                nadjc_obc(i) = nadjc_obc(i) + 1
                adjc_obc(i,nadjc_obc(i)) = ic
                IF(mod(type_obc(i),2) == 1)THEN
                   !!Node is Linear, Mark Cell as Linear for Flux Update
                   iucp(ic) = 0
                END IF
             END IF
          END DO
       END IF
    END DO

    DEALLOCATE(nxobc_tmp,nyobc_tmp)

    !--Determine Adjacent FluxFractions--------------------------------------------!

    ALLOCATE(nfluxf_obc(iobcn))  ; nfluxf_obc = 0
    ALLOCATE(fluxf_obc(iobcn,2)) ; fluxf_obc  = 0
    DO i=1,iobcn
       IF(nadjn_obc(i) == 1)THEN
          nfluxf_obc(i) = 1
          fluxf_obc(i,1) = 1.
          fluxf_obc(i,2) = 0.
       ELSE
          nfluxf_obc(i) = 2
          n1   = i_obc_n(i)
          n2   = adjn_obc(i,1)
          n3   = adjn_obc(i,2)
          e1 = sqrt( (vx(n1)-vx(n2))**2 + (vy(n1)-vy(n2))**2)
          e2 = sqrt( (vx(n1)-vx(n3))**2 + (vy(n1)-vy(n3))**2)
          fluxf_obc(i,1) =  e1/(e1+e2)
          fluxf_obc(i,2) =  e2/(e1+e2)
       END IF
    END DO
    !--Determine 1st Layer Neighbor for Open Boundary Points-----------------------!
    !--Node Chosen is Node That is Connected to OBC Node and is Oriented           !
    !--Most Normal to the Boundary.  It is not Necessarily the Closest Node        !
    !--Determine also DLTN_OBC, the normal component of the distance between       !
    !--Next_obc and the open boundary node                                         !

    DO i=1,iobcn
       dotmax =  -1.0
       inode = i_obc_n(i)
       DO j=1,ntsn(inode)-1
          jnode = nbsn(inode,j)
          IF(isonb(jnode) /= 2 .AND. inode /= jnode)THEN
             dxn_tmp = vx(jnode)-vx(inode)
             dyn_tmp = vy(jnode)-vy(inode)
             dxn = dxn_tmp/sqrt(dxn_tmp**2 + dyn_tmp**2)
             dyn = dyn_tmp/sqrt(dxn_tmp**2 + dyn_tmp**2)
             dot = dxn*nxobc(i) + dyn*nyobc(i)
             IF(dot > dotmax)THEN
                dotmax = dot
                next_obc(i) = jnode
             END IF
          END IF
       END DO
    END DO

    !--Determine 2nd Layer Neighbor for Open Boundary Points-----------------------!

    DO i=1,iobcn
       dotmax =  -1.0
       inode = next_obc(i)
       IF(inode > m) cycle
       DO j=1,ntsn(inode)-1
          jnode = nbsn(inode,j)
          IF(isonb(jnode) /= 2)THEN
             dxn_tmp = vx(jnode)-vx(inode)
             dyn_tmp = vy(jnode)-vy(inode)
             dxn = dxn_tmp/(sqrt(dxn_tmp**2 + dyn_tmp**2) + 1.0e-9_sp)
             dyn = dyn_tmp/(sqrt(dxn_tmp**2 + dyn_tmp**2) + 1.0e-9_sp)
             dot = dxn*nxobc(i) + dyn*nyobc(i)
             IF(dot > dotmax)THEN
                dotmax = dot
                next_obc2(i) = jnode
             END IF
          END IF
       END DO
    END DO

    !--Determine DLTN_OBC----------------------------------------------------------!
    ALLOCATE(dltn_obc(iobcn))
    DO i=1,iobcn
       i1 = i_obc_n(i)
       i2 = next_obc(i)

       dx = vx(i2)-vx(i1)
       dy = vy(i2)-vy(i1)
       dltn_obc(i) = abs(dx*nxobc(i) + dy*nyobc(i))
    END DO

!!$   RETURN
!!$!--Dump Information to Matlab Files for Checking-------------------------------!
!!$
!!$   OPEN(UNIT=81,FILE='mesh.scatter',FORM='formatted')
!!$   DO I=1,M
!!$     WRITE(81,*)vx(i),vy(i)
!!$   END DO
!!$   CLOSE(81)
!!$   OPEN(UNIT=81,FILE='nextobc.scatter',FORM='formatted')
!!$   DO I=1,IOBCN
!!$     I1 = NEXT_OBC(I)
!!$     WRITE(81,*)VX(I1),VY(I1)
!!$   END DO
!!$   CLOSE(81)
!!$   OPEN(UNIT=81,FILE='nextobc2.scatter',FORM='formatted')
!!$   DO I=1,IOBCN
!!$     I1 = NEXT_OBC2(I)
!!$     WRITE(81,*)VX(I1),VY(I1)
!!$   END DO
!!$   CLOSE(81)
!!$   OPEN(UNIT=81,FILE='iobcn.scatter',FORM='formatted')
!!$   DO I=1,IOBCN
!!$     I1 = I_OBC_N(I)
!!$     WRITE(81,*)VX(I1),VY(I1) 
!!$   END DO
!!$   CLOSE(81)
!!$   OPEN(UNIT=81,FILE='obcnorm.scatter',FORM='formatted')
!!$   DO I=1,IOBCN
!!$     I1 = I_OBC_N(I)
!!$     WRITE(81,*)NXOBC(I1),NYOBC(I1) 
!!$   END DO
!!$   CLOSE(81)
!!$   OPEN(UNIT=81,FILE='nonlinear.scatter',FORM='formatted')
!!$   DO I=1,IOBCN
!!$     IF(NADJC_OBC(I) > 0) WRITE(81,*)XC(ADJC_OBC(I,1)),YC(ADJC_OBC(I,1)) 
!!$     IF(NADJC_OBC(I) > 1) WRITE(81,*)XC(ADJC_OBC(I,2)),YC(ADJC_OBC(I,2)) 
!!$   END DO
!!$   CLOSE(81)
!!$   OPEN(UNIT=81,FILE='linear.scatter',FORM='formatted')
!!$   DO I=1,N
!!$     IF(IUCP(I)==0)THEN
!!$       WRITE(81,*)XC(I),YC(I) 
!!$     END IF
!!$   END DO
!!$   CLOSE(81)

    RETURN
  END SUBROUTINE setup_obc
  !==============================================================================|

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

  !==============================================================================|
  SUBROUTINE flux_obn(K)

    USE all_vars
    IMPLICIT NONE

    INTEGER, INTENT(IN)  :: K
    INTEGER              :: I,J,C1,C2
    REAL(SP)             :: FF,FLUX 

    fluxobn = 0.0_sp

    DO i = 1, iobcn

       j = i_obc_n(i)
       !Compute Boundary Flux From Continuity Flux Defect
       flux = -(elf(j)-elrk(j))*art1(j)/(alpha_rk(k)*dte)-xflux_obcn(i)

       !Set Flux In Adjacent Nonlinear BC Element 1 (If Exists)
       IF(nadjc_obc(i) > 0) THEN
          c1 = adjc_obc(i,1) 
          ff = fluxf_obc(i,1)
          fluxobn(c1) = fluxobn(c1) + ff*flux 
       END IF

       !Set Flux In Adjacent Nonlinear BC Element 2 (If Exists)
       IF(nadjc_obc(i) > 1) THEN
          c2 =  adjc_obc(i,2) 
          ff = fluxf_obc(i,2)
          fluxobn(c2) = fluxobn(c2) + ff*flux 
       END IF

    END DO

    RETURN
  END SUBROUTINE flux_obn
  !========================================================================

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

  !==============================================================================|
  SUBROUTINE bcond_t_perturbation(T2D_NEXT,T2D,TTMP,I,J,J1)
    !==============================================================================|
    ! Calculate the OBC for temperature perturbation                               |
    !==============================================================================|
    USE all_vars
    IMPLICIT NONE

    !   INTEGER :: I1,I2,J,JN
    INTEGER :: I,J,J1,J2,K
    REAL(SP):: CC,CP,MU,CL
    REAL(SP):: PERT_NEXT,PERT,T2D_NEXT,T2D
    REAL(SP):: T2D_NEXT1,TM12D_NEXT2,TM12D_NEXT1,TM22D_NEXT1
    REAL(SP):: TTMP(IOBCN,KBM1)

    SELECT CASE(type_tsobc)

    CASE(1)
       DO k=1,kbm1
          ttmp(i,k) = tf1(j1,k) - t2d_next
       END DO
    CASE(2)
       cc = sqrt(grav_n(j)*h(j))*dti/dltn_obc(i)
       cp = cc + 1.0_sp
       DO k=1,kbm1
          pert_next = tf1(j1,k) - t2d_next
          pert      = t1(j,k) - t2d
          ttmp(i,k) = (cc*pert_next + pert)/cp
       END DO
    CASE(3)
       cc = sqrt(grav_n(j)*h(j))*dti/dltn_obc(i)
       cp = cc + 1.0_sp
       DO k=1,kbm1
          pert_next = tf1(j1,k) - t2d_next
          pert      = t1(j,k) - t2d
          ttmp(i,k) = (cc*pert_next + pert*(1.0_sp - dti/10800.0_sp))/cp
       END DO
    CASE(4)
       j2 = next_obc2(i)
       t2d_next1  =0.0_sp
       tm12d_next2=0.0_sp
       tm12d_next1=0.0_sp
       tm22d_next1=0.0_sp
       DO k=1,kbm1
          t2d_next1  =t2d_next1  +t1(j1,k)*dz(j1,k)
          tm12d_next2=tm12d_next2+t1m1(j2,k)*dz(j2,k)
          tm12d_next1=tm12d_next1+t1m1(j,k)*dz(j,k)
          tm22d_next1=tm22d_next1+t1m2(j1,k)*dz(j1,k)
       END DO

       DO k=1,kbm1
          cl = ((t1m2(j1,k)-tm22d_next1)-(t1(j1,k)-t2d_next1))/   &
               ((t1(j1,k)-t2d_next1)+(t1m2(j1,k)-tm22d_next1)     &
               -2.0*(t1m1(j2,k)-tm12d_next2))
          IF(cl >= 1.0)THEN
             mu = 1.0
          ELSE IF(cl > 0.0 .AND. cl < 1.0)THEN
             mu = cl
          ELSE
             mu = 0.0
          END IF

          ttmp(i,k)=((t1m1(j,k)-tm12d_next1)*(1.0-mu)     &
               +2.0*mu*(t1(j1,k)-t2d_next1))/(1.0+mu)
       END DO

    CASE DEFAULT
       CALL fatal_error("INVALID OBC_TS_TYPE IN NML_OPEN_BOUNDARY_CONTROL"&
            &, "See mod_obcs.F")

    END SELECT

    RETURN
  END SUBROUTINE bcond_t_perturbation
  !========================================================================

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

  !==============================================================================|
  SUBROUTINE bcond_s_perturbation(S2D_NEXT,S2D,STMP,I,J,J1)
    !==============================================================================|
    ! Calculate the OBC for salinity perturbation                                  |
    !==============================================================================|
    USE all_vars
    IMPLICIT NONE

    INTEGER :: I,J,J1,J2,K
    REAL(SP):: CC,CP,MU,CL
    REAL(SP):: PERT_NEXT,PERT,S2D_NEXT,S2D
    REAL(SP):: S2D_NEXT1,SM12D_NEXT2,SM12D_NEXT1,SM22D_NEXT1
    REAL(SP):: STMP(IOBCN,KBM1)

    SELECT CASE(type_tsobc)

    CASE(1)
       DO k=1,kbm1
          stmp(i,k) = sf1(j1,k) - s2d_next
       END DO
    CASE(2)
       cc = sqrt(grav_n(j)*h(j))*dti/dltn_obc(i)
       cp = cc + 1.0_sp
       DO k=1,kbm1
          pert_next = sf1(j1,k) - s2d_next
          pert      = s1(j,k) - s2d
          stmp(i,k) = (cc*pert_next + pert)/cp
       END DO
    CASE(3)
       cc = sqrt(grav_n(j)*h(j))*dti/dltn_obc(i)
       cp = cc + 1.0_sp
       DO k=1,kbm1
          pert_next = sf1(j1,k) - s2d_next
          pert      = s1(j,k) - s2d
          stmp(i,k) = (cc*pert_next + pert*(1.0_sp - dti/10800.0_sp))/cp
       END DO
    CASE(4)
       j2 = next_obc2(i)
       s2d_next1  =0.0_sp
       sm12d_next2=0.0_sp
       sm12d_next1=0.0_sp
       sm22d_next1=0.0_sp
       DO k=1,kbm1
          s2d_next1  =s2d_next1  +s1(j1,k)*dz(j1,k)
          sm12d_next2=sm12d_next2+s1m1(j2,k)*dz(j2,k)
          sm12d_next1=sm12d_next1+s1m1(j,k)*dz(j,k)
          sm22d_next1=sm22d_next1+s1m2(j1,k)*dz(j1,k)
       END DO

       DO k=1,kbm1
          cl = ((s1m2(j1,k)-sm22d_next1)-(s1(j1,k)-s2d_next1))/   &
               ((s1(j1,k)-s2d_next1)+(s1m2(j1,k)-sm22d_next1)     &
               -2.0*(s1m1(j2,k)-sm12d_next2))
          IF(cl >= 1.0)THEN
             mu = 1.0
          ELSE IF(cl > 0.0 .AND. cl < 1.0)THEN
             mu = cl
          ELSE
             mu = 0.0
          END IF

          stmp(i,k)=((s1m1(j,k)-sm12d_next1)*(1.0-mu)     &
               +2.0*mu*(s1(j1,k)-s2d_next1))/(1.0+mu)
       END DO

    CASE DEFAULT
       CALL fatal_error("INVALID OBC_TS_TYPE IN NML_OPEN_BOUNDARY_CONTROL"&
            &, "See mod_obcs.F")

    END SELECT

    RETURN
  END SUBROUTINE bcond_s_perturbation
  !========================================================================

  SUBROUTINE gday1(IDD,IMM,IYY,ICC,KD)
!
!  given day,month,year and century(each 2 digits), gday returns
!  the day#, kd based on the gregorian calendar.
!  the gregorian calendar, currently 'universally' in use was
!  initiated in europe in the sixteenth century. note that gday
!  is valid only for gregorian calendar dates.
!
!  kd=1 corresponds to january 1, 0000
!
!  note that the gregorian reform of the julian calendar 
!  omitted 10 days in 1582 in order to restore the date
!  of the vernal equinox to march 21 (the day after
!  oct 4, 1582 became oct 15, 1582), and revised the leap 
!  year rule so that centurial years not divisible by 400
!  were not leap years.
!
!  this routine was written by eugene neufeld, at ios, in june 1990.
!
    integer idd, imm, iyy, icc, kd
    integer ndp(13)
    integer ndm(12)
    data ndp/0,31,59,90,120,151,181,212,243,273,304,334,365/
    data ndm/31,28,31,30,31,30,31,31,30,31,30,31/
!
!  test for invalid input:
    if(icc.lt.0)then
!     write(11,5000)icc
      call pstop
    endif
    if(iyy.lt.0.or.iyy.gt.99)then
!     write(11,5010)iyy
      call pstop
    endif
    if(imm.le.0.or.imm.gt.12)then
!     write(11,5020)imm
      call pstop
      endif
    if(idd.le.0)then
!     write(11,5030)idd
      call pstop
    endif
    if(imm.ne.2.and.idd.gt.ndm(imm))then
!     write(11,5030)idd
      call pstop
    endif
    if(imm.eq.2.and.idd.gt.29)then
!     write(11,5030)idd
      call pstop
    endif
    if(imm.eq.2.and.idd.gt.28.and.((iyy/4)*4-iyy.ne.0.or.(iyy.eq.0.and.(icc/4)*4-icc.ne.0)))then
!     write(11,5030)idd
      call pstop
    endif
5000  format(' input error. icc = ',i7)
5010  format(' input error. iyy = ',i7)
5020  format(' input error. imm = ',i7)
5030  format(' input error. idd = ',i7)
!
!  calculate day# of last day of last century:
    kd = icc*36524 + (icc+3)/4
!
!  calculate day# of last day of last year:
    kd = kd + iyy*365 + (iyy+3)/4
!
!  adjust for century rule:
!  (viz. no leap-years on centurys except when the 2-digit
!  century is divisible by 4.)
    if(iyy.gt.0.and.(icc-(icc/4)*4).ne.0) kd=kd-1
!  kd now truly represents the day# of the last day of last year.
!
!  calculate day# of last day of last month:
    kd = kd + ndp(imm)
!
!  adjust for leap years:
    if(imm.gt.2.and.((iyy/4)*4-iyy).eq.0.and.((iyy.ne.0).or.(((icc/4)*4-icc).eq.0))) kd=kd+1
!  kd now truly represents the day# of the last day of the last
!  month.
!
!  calculate the current day#:
    kd = kd + idd

  RETURN
  END SUBROUTINE gday1 

  subroutine astro(d1,h,pp,s,p,np,dh,dpp,ds,dp,dnp)
! this subroutine calculates the following five ephermides
! of the sun and moon
! h = mean longitude of the sum
! pp = mean longitude of the solar perigee
! s = mean longitude of the moon
! p = mean longitude of the lunar perigee
! np = negative of the longitude of the mean ascending node
! and their rates of change.
! units for the ephermides are cycles and for their derivatives
! are cycles/365 days
! the formulae for calculating this ephermides were taken from
! pages 98 and 107 of the explanatory supplement to the
! astronomical ephermeris and the american ephermis and
! nautical almanac (1961)
!
  implicit real*8(a-h,o-z)
  real*8 np
  d2=d1*1.d-4
  f=360.d0
  f2=f/365.d0
  h=279.696678d0+.9856473354d0*d1+.00002267d0*d2*d2
  pp=281.220833d0+.0000470684d0*d1+.0000339d0*d2*d2+.00000007d0*d2**3
  s=270.434164d0+13.1763965268d0*d1-.000085d0*d2*d2+.000000039d0*d2**3
  p=334.329556d0+.1114040803d0*d1-.0007739d0*d2*d2-.00000026d0*d2**3
  np=-259.183275d0+.0529539222d0*d1-.0001557d0*d2*d2-.00000005d0*d2**3
  h=h/f
  pp=pp/f
  s=s/f
  p=p/f
  np=np/f
  h=h-dint(h)
  pp=pp-dint(pp)
  s=s-dint(s)
  p=p-dint(p)
  np=np-dint(np)
  dh=.9856473354d0+2.d-8*.00002267d0*d1
  dpp=.0000470684d0+2.d-8*.0000339d0*d1+3.d-12*.00000007d0*d1**2
  ds=13.1763965268d0-2.d-8*.000085d0*d1+3.d-12*.000000039d0*d1**2
  dp=.1114040803d0-2.d-8*.0007739d0*d1-3.d-12*.00000026d0*d1**2
  dnp=+.0529539222d0-2.d-8*.0001557d0*d1-3.d-12*.00000005d0*d1**2
  dh=dh/f2
  dpp=dpp/f2
  ds=ds/f2
  dp=dp/f2
  dnp=dnp/f2
  return
  end subroutine astro

!!$

END MODULE mod_obcs