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


!==============================================================================|
!     SET BOUNDARY CONDITIONS FOR ALMOST ALL VARIABLES                         |
!                                                                              |
!         idx: identifies which variables are considered                       |
!              1=tidal forcing                                                 |
!              2=solid bcs for external mode uaf and vaf                       |
!              3=solid bcs for internal mode uf and vf                         |
!              4=open bcs for s and t                                          |
!              5=solid bcs for internal mode u and v                           |
!              6=unused                                                        |
!              7=unused                                                        |
!              8=the surface forcings for internal mode                        |
!              9=the surface forcings for external mode                        |
!                                                                              |
!==============================================================================|

   SUBROUTINE bcond_gcn(IDX,K_RK)

!==============================================================================|
   USE all_vars
   USE bcs
   USE mod_obcs
   USE mod_force
   USE mod_par
   USE mod_wd
   USE mod_bulk




   IMPLICIT NONE
   INTEGER, INTENT(IN) :: IDX
   REAL(SP) :: ZREF(KBM1),ZREFJ(KBM1),TSIGMA(KBM1),SSIGMA(KBM1)
   REAL(SP) :: TTMP(KBM1),STMP(KBM1),TREF(KBM1),SREF(KBM1)
   REAL(SP) :: PHY_Z(KBM1),PHY_Z1(KBM1)
   REAL(SP) :: TT1(KBM1),TT2(KBM1),SS1(KBM1),SS2(KBM1)
!   REAL(SP) :: TIME1,FACT,UFACT,FORCE,QPREC,QEVAP,UI,VI,UNTMP,VNTMP,TX,TY,HFLUX
   REAL(SP) :: TIME1,FACT,UFACT,FORCE,UI,VI,UNTMP,VNTMP,TX,TY,HFLUX
!   REAL(SP) :: DTXTMP,DTYTMP,QPREC2,QEVAP2,SPRO,WDS,CD,SPCP,ROSEA
   REAL(SP) :: DTXTMP,DTYTMP,SPRO,WDS,CD,SPCP,ROSEA,ROSEA1(MT),SPRO1(MT)
   REAL(SP) :: PHAI_IJ,ALPHA1,DHFLUXTMP,DHSHORTTMP,HSHORT,TIMERK1

   ! VARIABLES FOR LONG SHORE FLOW STUFF
   REAL(SP) :: ANG_WND,WNDALONG,RHOINTNXT,RHOINTCUR,CUMEL
   REAL(SP) :: TAU_X,TAU_Y, mag_wnd
   REAL(SP), POINTER,DIMENSION(:) :: lcl_dat, gbl_dat


   REAL(SP),POINTER :: eta_lcl(:), eta_gbl(:) ,elfgeo_gbl(:),elfgeo_lcl(:)

   INTEGER  I,J,K,I1,I2,J1,J2,II,L1,L2,IERR
   INTEGER  cdx,ndx,KDAM_TMP,K_RK

!!$!=================|DEBUG|============================
!!$   LOGICAL, save :: INIT=.false.
!!$   INTEGER, save :: icount
!!$   character(len=5):: ccount
!!$   REAL(SP),POINTER :: temp_lcl(:),temp_gl(:) 
!!$!=================|DEBUG|============================


  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: bcond_gcn: ",idx


   SELECT CASE(idx)
!==============================================================================|
   CASE(1) !Surface Elevation Boundary Conditions (Tidal Forcing)              !
!==============================================================================|
   IF(obc_elevation_forcing_on) CALL bcond_asl
   CALL bcond_asl_clp
   CALL bcond_gwi(k_rk)
   CALL bcond_bki(k_rk)
   CALL bcond_ore

!
!--Allow setup/down on north boundary in response to longshore wind
!--Corrects for Wind-Driven Barotropic Response (See Schwing 1989)
!--Implemented by Jamie Pringle - Rewritten for parallel by D.Stuebe
!

   IF (obc_longshore_flow_on) THEN
      
      ! CALCULATE THE ADJUSTMENT DUE TO WIND DRIVEN FLOW ALLONG THE
      ! COAST OF NOVA SCOTIA
      DO i = 1,nobclsf
         tau_x = 0.0_sp
         tau_y = 0.0_sp

         i1= ibclsf(i)
         DO j = 1,ntve(i1)
            k = nbve(i1,j)
            tau_x = tau_x + wusurf2(k)
            tau_y = tau_y + wvsurf2(k)
         END DO
         
         tau_x = tau_x/real(ntve(i1),sp)
         tau_y = tau_y/real(ntve(i1),sp)

         tau_x = tau_x * 1000.0_sp ! Approximate scale factor
         tau_y = tau_y * 1000.0_sp ! Approximate scale factor

         mag_wnd = sqrt(tau_x**2+tau_y**2)
         
         
         IF (mag_wnd .GT. 0.0_sp) THEN
            ang_wnd=atan2(tau_y,tau_x) - wdf_ang(i)
         ELSE
            ang_wnd=0.0_sp
         END IF
         

         wndalong=sin(ang_wnd)*mag_wnd

         ! SUBTRACT THE COMPONET ALONG SHORE AND INTO THE DOMAIN FROM
         ! THE SEA SURFACE HEIGHT
         
         elf(i1)=elf(i1) - wndalong*rbc_wdf(i)

!!$         !=================|DEBUG|============================
!!$ THIS IS NOT PARALLEL SAFE!
!!$            if(.not. INIT)then
!!$               init =.true.
!!$               icount=0
!!$            end if
!!$            
!!$            IF(I==1) THEN
!!$               write(ccount,'(I5.5)') icount
!!$               
!!$               call FOPEN(157,"WDF_"//ccount,'ofr')
!!$               write(157,*) "WDF DEBUG: ",icount
!!$               write(157,*) "Node#, ANG,  ELFGEO"
!!$            END IF
!!$
!!$            write(157,*) IBCLSF(I), WDF_ANG(I)/deg2rad, RBC_WDF(I),&
!!$                 &-1.0_SP*WNDALONG*RBC_WDF(I), MAG_WND*SIN(ANG_WND)
!!$
!!$            IF(I==nobclsf) THEN
!!$               icount = icount +1
!!$               
!!$               close(157)
!!$            END IF
!!$!=================|DEBUG|============================

      END DO
      
      

      ! NOW ADJUSTMENT ELF UNDER THERMAL WIND SUCH THAT THE BOTTOM
      ! VELOCITY IS ZERO 
      allocate(eta_lcl(nobclsf)); eta_lcl=0.0_sp
! DEBUG
!      allocate(temp_lcl(NOBCLSF)); temp_lcl=0.0_sp

      rhointcur= 0.0_sp
      rhointnxt = 0.0_sp

      DO i=nobclsf,1,-1  ! COUNT BACKWARDS - The onshore direction
         ndx  = nbclsf(i)
         cdx = ibclsf(i)

         !INTEGRATE RHO IN DEPTH, CONVERT TO MKS UNITS DAMIT
         rhointcur=rhointnxt
         rhointnxt=0.0_sp
         DO k=1,kbm1
            rhointnxt=rhointnxt+(1.0_sp+rho1(cdx,k))*1.0e3_sp*dz(cdx,k)
         END DO
! DEBUG         
!         temp_lcl(I)=RHOINTNXT

!ESTIMATE DENSITY GRADIENT, AND MODIFY BOUNDARY ELEVATION
!NOTE THE FACTOR OF 1000 AND 2 TO COMPENSATE FOR THE
!FACT THAT THE MODEL STORES RHO1 AS SIGMA, AND IN CGS.

! ADJUST THE SEA SURFACE HEIGHT FOR A MEAN FLOW DUE TO DENSITY GRADIENT
! CALCULATE THE THERMAL WIND AND ADJUST THE SEA SURFACE HEIGHT SO THAT THE BOTTOM
! BOUNDARY IS A LEVEL OF NO MOTION:
!
! ETATAN=[ -1/(Rho_bar)*(del(h*rho)/del(s) - Rho_bot*(del(h)/del(s) ] *del(s)
!
         IF (i /= nobclsf) THEN
            eta_lcl(i)=-(1.0_sp/(0.5_sp*(rhointnxt+rhointcur))) &
                 *((h(cdx)*rhointnxt-h(ndx)*rhointcur) &
                 -0.5_sp*1.0e3_sp*(2.0_sp+rho1(cdx,kbm1)+rho1(ndx,kbm1)) &
                 *(h(cdx)-h(ndx)))
         END IF

      END DO

      cumel = 0.0_sp
      
      IF(serial) THEN
         DO i=nobclsf,1,-1 ! COUNT BACKWARD FROM THE SHELF TOWARD SHORE
            ndx=ibclsf(i)
            cumel=cumel+eta_lcl(i)
            elf(ndx)=elf(ndx)+ cumel *rbc_geo(i)
!            write(ipt,*) "NDX=",NDX,"ELF(NDX)=",ELF(NDX)
         END DO
      END IF

! DEBUG
!      deallocate(temp_lcl)
      DEALLOCATE(eta_lcl)

   END IF


!==============================================================================|
   CASE(2) !External Mode Velocity Boundary Conditions                         |
!==============================================================================|

   DO i=1,n

!
!--2 SOLID BOUNDARY EDGES------------------------------------------------------|
!
   IF(isbce(i) == 3) THEN
     uaf(i)=0.0_sp
     vaf(i)=0.0_sp
   END IF

!  GWC SPEED REPLACE ABOVE 4 LINES
!   UAF(LISBCE_3(1:NISBCE_3)) = 0.  
!   VAF(LISBCE_3(1:NISBCE_3)) = 0.  

!
!--1 SOLID BOUNDARY EDGE-------------------------------------------------------|
!
   IF(isbce(i) == 1) THEN
     alpha1=alpha(i)
     IF(numqbc > 0) THEN
       IF(river_inflow_location == 'node') THEN
         DO j=1,numqbc
           i1=inodeq(j)
           j1=nbve(i1,1)
           j2=nbve(i1,ntve(i1))
           IF((i == j1).OR.(i == j2)) THEN
             untmp=uaf(i)*cos(angleq(j))+vaf(i)*sin(angleq(j))
             vntmp=-uaf(i)*sin(angleq(j))+vaf(i)*cos(angleq(j))
             untmp=max(untmp,0.0_sp)
             uaf(i)=untmp*cos(angleq(j))-vntmp*sin(angleq(j))
             vaf(i)=untmp*sin(angleq(j))+vntmp*cos(angleq(j))
               GOTO 21
               END IF 
             END DO
            ELSE IF(river_inflow_location == 'edge') THEN
             DO j=1,numqbc
               j1=icellq(j)
               IF(i == j1) THEN
                untmp=uaf(i)*cos(angleq(j))+vaf(i)*sin(angleq(j))
                vntmp=-uaf(i)*sin(angleq(j))+vaf(i)*cos(angleq(j))
                untmp=max(untmp,0.0_sp)
                uaf(i)=untmp*cos(angleq(j))-vntmp*sin(angleq(j))
                vaf(i)=untmp*sin(angleq(j))+vntmp*cos(angleq(j))
               GOTO 21
               END IF
             END DO
           END IF
            
           END IF 


           ui= uaf(i)*(sin(alpha1))**2-vaf(i)*sin(alpha1)*cos(alpha1)
           vi=-uaf(i)*sin(alpha1)*cos(alpha1)+vaf(i)*(cos(alpha1))**2
           uaf(i)=ui
           vaf(i)=vi

!!#          if defined (1)
!!           UI= UAS(I)*(SIN(ALPHA1))**2-VAS(I)*SIN(ALPHA1)*COS(ALPHA1)
!!           VI=-UAS(I)*SIN(ALPHA1)*COS(ALPHA1)+VAS(I)*(COS(ALPHA1))**2
!!           UAS(I)=UI
!!           VAS(I)=VI
!!#          endif


21        CONTINUE
          END IF
       END DO

!!# if defined (THIN_DAM)
!!    DO I=1,NT
!!      IF(CLP_CELL(I)==1)THEN
!!        ALPHA1=CLP_ALPHA(I)
!!        UNTMP=-UAF(I)*COS(ALPHA1)-VAF(I)*SIN(ALPHA1)
!!        VNTMP=UAF(I)*SIN(ALPHA1)-VAF(I)*COS(ALPHA1)
!!        IF(CLP_KDAM(I,1)==1)KDAM_TMP=KDAM1(CLP_KDAM(I,2))
!!        IF(CLP_KDAM(I,1)==2)KDAM_TMP=(KDAM1(CLP_KDAM(I,2))+KDAM1(CLP_KDAM(I,2)))/2
!!        IF(KDAM_TMP<=KBM1/2)UNTMP=UNTMP*KDAM_TMP/KBM1
!!        UAF(I)=-UNTMP*COS(ALPHA1)+VNTMP*SIN(ALPHA1)
!!        VAF(I)=-UNTMP*SIN(ALPHA1)-VNTMP*COS(ALPHA1)
!!      END IF
!!    END DO
!!# endif

!==============================================================================|
   CASE(3) !3-D Velocity Boundary Conditions                                   !
!==============================================================================|

!  GWC SPEED REPLACE NEXT 4 LINES
!   UF(LISBCE_3(1:NISBCE_3),1:KBM1) = 0.  
!   VF(LISBCE_3(1:NISBCE_3),1:KBM1) = 0.  

      ! Suggest move the kbm1 loop inside 'if((i.eq.j1).or.(i.eq.j2)) then'

      ! THIS SHOULD BE SAFE FOR MORE THAN ONE RIVER AT THE SAME NODE
      ! OR EDGE

       do i= 1, n
         do k =1, kbm1
          if(isbce(i).eq.3) then
           uf(i,k)=0.0_sp
           vf(i,k)=0.0_sp
          end if

          if(isbce(i).eq.1) then
          alpha1=alpha(i)
          if(numqbc.ge.1) then
          if(river_inflow_location.eq.'node') then
            do j=1,numqbc
              i1=inodeq(j)
              j1=nbve(i1,1)
              j2=nbve(i1,ntve(i1))
              if((i.eq.j1).or.(i.eq.j2)) then
               untmp=uf(i,k)*cos(angleq(j))+vf(i,k)*sin(angleq(j))
               vntmp=-uf(i,k)*sin(angleq(j))+vf(i,k)*cos(angleq(j))
               untmp=max(untmp,0.0_sp)
               uf(i,k)=untmp*cos(angleq(j))-vntmp*sin(angleq(j))
               vf(i,k)=untmp*sin(angleq(j))+vntmp*cos(angleq(j))
              goto 31
              end if
            end do
           else if(river_inflow_location.eq.'edge') then
            do j=1,numqbc
              j1=icellq(j)
              if(i.eq.j1) then
               untmp=uf(i,k)*cos(angleq(j))+vf(i,k)*sin(angleq(j))
               vntmp=-uf(i,k)*sin(angleq(j))+vf(i,k)*cos(angleq(j))
               untmp=max(untmp,0.0_sp)
               uf(i,k)=untmp*cos(angleq(j))-vntmp*sin(angleq(j))
               vf(i,k)=untmp*sin(angleq(j))+vntmp*cos(angleq(j))
              goto 31
              end if
            end do
           else
             print*, 'river_inflow_location not correct'
             call pstop
          end if
          end if


          ui= uf(i,k)*(sin(alpha1))**2-vf(i,k)*sin(alpha1)*cos(alpha1)
          vi=-uf(i,k)*sin(alpha1)*cos(alpha1)+vf(i,k)*(cos(alpha1))**2
          uf(i,k)=ui
          vf(i,k)=vi

!!#         if defined (1)
!!          ui= us(i,k)*(sin(alpha1))**2-vs(i,k)*sin(alpha1)*cos(alpha1)
!!          vi=-us(i,k)*sin(alpha1)*cos(alpha1)+vs(i,k)*(cos(alpha1))**2
!!          us(i,k)=ui
!!          vs(i,k)=vi
!!#         endif


31        continue
          end if
         end do
       end do



!==============================================================================|
   CASE(4) !Blank                                                              !
!==============================================================================|

!==============================================================================|
   CASE(5) !!SOLID BOUNDARY CONDITIONS ON U AND V                              !
!==============================================================================|


!  GWC SPEED REPLACE NEXT 4 LINES
!   U(LISBCE_3(1:NISBCE_3),1:KBM1) = 0.  
!   V(LISBCE_3(1:NISBCE_3),1:KBM1) = 0.  

       do i= 1, n
         do k =1, kbm1
          if(isbce(i).eq.3) then
          u(i,k)=0.0_sp
          v(i,k)=0.0_sp
          end if

          if(isbce(i).eq.1) then
          alpha1=alpha(i)
          if(numqbc.ge.1) then
          if(river_inflow_location.eq.'node') then
            do j=1,numqbc
              i1=inodeq(j)
              j1=nbve(i1,1)
              j2=nbve(i1,ntve(i1))
              if((i.eq.j1).or.(i.eq.j2)) then
               untmp=u(i,k)*cos(angleq(j))+v(i,k)*sin(angleq(j))
               vntmp=-u(i,k)*sin(angleq(j))+v(i,k)*cos(angleq(j))
               untmp=max(untmp,0.0_sp)
               u(i,k)=untmp*cos(angleq(j))-vntmp*sin(angleq(j))
               v(i,k)=untmp*sin(angleq(j))+vntmp*cos(angleq(j))
              goto 51
              end if
            end do
           else if(river_inflow_location.eq.'edge') then
            do j=1,numqbc
              j1=icellq(j)
              if(i.eq.j1) then
               untmp=u(i,k)*cos(angleq(j))+v(i,k)*sin(angleq(j))
               vntmp=-u(i,k)*sin(angleq(j))+v(i,k)*cos(angleq(j))
               untmp=max(untmp,0.0_sp)
               u(i,k)=untmp*cos(angleq(j))-vntmp*sin(angleq(j))
               v(i,k)=untmp*sin(angleq(j))+vntmp*cos(angleq(j))
              goto 51
              end if
            end do
           else
             print*, 'river_inflow_location not correct'
             call pstop
          end if
          end if


          ui= u(i,k)*(sin(alpha1))**2-v(i,k)*sin(alpha1)*cos(alpha1)
          vi=-u(i,k)*sin(alpha1)*cos(alpha1)+v(i,k)*(cos(alpha1))**2
          u(i,k)=ui
          v(i,k)=vi

!!#         if defined (1)
!!          ui= us(i,k)*(sin(alpha1))**2-vs(i,k)*sin(alpha1)*cos(alpha1)
!!          vi=-us(i,k)*sin(alpha1)*cos(alpha1)+vs(i,k)*(cos(alpha1))**2
!!          us(i,k)=ui
!!          vs(i,k)=vi
!!#         endif


51        continue
          end if
         end do
       end do



!==============================================================================|
   CASE(6) !Blank                                                              !
!==============================================================================|

!==============================================================================|
   CASE(7) !Blank                                                              !
!==============================================================================|

!==============================================================================|
   CASE(8) !!SURFACE FORCING FOR INTERNAL MODE                                 !
!==============================================================================|

!
!--Fresh Water Discharge-------------------------------------------------------|
!

      IF(numqbc_gl .GT. 0) THEN
         CALL update_rivers(inttime,qdis,temp=tdis,salt=sdis)
                  
         qdis    = qdis*ramp
         
!#  if defined (WATER_QUALITY)
!         DO N1 = 1, NB
!            ! NOT YET OPERATIONAL!!!
!            WDIS(:,N1) = UFACT*DWDIS(:,N1,L1) + FACT*DWDIS(:,N1,L2)
!         END DO
!#  endif
         
      END IF
      
      
      
      IF (wind_on) THEN
         IF (wind_type == speed)THEN
            CALL update_wind(inttime,uuwind,vvwind)
            
            CALL asimple_drag(uuwind,vvwind,wusurf,wvsurf)

         ELSEIF(wind_type == stress)THEN
            CALL update_wind(inttime,wusurf,wvsurf)
         END IF
         

         ! For output only - don't divide by density
         wusurf_save = wusurf * ramp
         wvsurf_save = wvsurf * ramp

         ! Divide by density 
         wusurf = -wusurf * ramp *0.001_sp
         wvsurf = -wvsurf * ramp *0.001_sp
         
         ! MAJOR MISTAKE: NEED THE NEGATIVE SIGN FOR THE INTERNAL
         ! STEP WIND STRESS

      END IF
      
      IF (heating_on) THEN
         CALL update_heat(inttime,swrad_watts,wtsurf_watts)
         
         ! LOAD INTO THES VARIABLE TO SAVE THE FORCING IN THE CORRECT UNITS
         swrad_watts  = swrad_watts * ramp
         wtsurf_watts = wtsurf_watts * ramp

         spcp=4.2174e3_sp
         rosea = 1.023e3_sp
         spro = spcp*rosea
         
         wtsurf    = -wtsurf_watts/spro
         swrad     = -swrad_watts/spro


      END IF
      
      
      IF (precipitation_on) THEN
         CALL update_precipitation(inttime,qprec,qevap)
         ! NO RAMP FOR PRECIP/EVAP
      END IF






!
!-- Set Groundwater flux ------------------------------------------------------|
!
      
      IF (groundwater_on) THEN
         CALL update_groundwater(inttime,bfwdis,gw_temp=bfwtmp,gw_salt=bfwslt)
         bfwdis = ramp * bfwdis
         ! DO NOT RAMP TEMP AND SALINITY
      END IF
      
!==============================================================================|
   CASE(9) !External Mode Surface BCs (River Flux/Wind Stress/Heat/Moist)      !
!==============================================================================|
!
!-Freshwater Flux: Set  Based on Linear Interpolation Between Two Data Times---|
!

      IF(numqbc_gl .GT. 0) THEN
         CALL update_rivers(exttime,qdis2)
         
         qdis2 = qdis2*ramp
         
      END IF
      
!
!-- Set Precipitation/Evaporation/Surface Wind ---------------------------------|
!

      IF (precipitation_on) THEN
         CALL update_precipitation(exttime,qprec2,qevap2)
         ! NO RAMP FOR PRECIP/EVAP
      END IF
      
      IF (wind_on) THEN

         IF (wind_type == speed)THEN
            CALL update_wind(exttime,uuwind,vvwind)
            
            CALL asimple_drag(uuwind,vvwind,wusurf2,wvsurf2)

         ELSEIF(wind_type == stress)THEN
            CALL update_wind(exttime,wusurf2,wvsurf2)
         END IF


         wusurf2 = wusurf2 * ramp *0.001_sp
         wvsurf2 = wvsurf2 * ramp *0.001_sp
         
      END IF
   
!
!-- Set Groundwater flux ------------------------------------------------------|
!

      IF (groundwater_on) THEN
         CALL update_groundwater(exttime,bfwdis2)
         bfwdis2 = ramp * bfwdis2
      END IF
   END SELECT
   
   if(dbg_set(dbg_sbr)) write(ipt,*)&
        & "End: bcond_gcn"
   
 END SUBROUTINE bcond_gcn
!==============================================================================|