bcond_gcy.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	bcond_gcy (IDX, K_RK)

Function/Subroutine Documentation

bcond_gcy()

subroutine bcond_gcy	(	integer, intent(in)	IDX,
		integer	K_RK
	)

Definition at line 57 of file bcond_gcy.f90.

!==============================================================================|
   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,K_RK


  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: bcond_gcy: ",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)

      END DO
      
      

      ! NOW ADJUSTMENT ELF UNDER THERMAL WIND SUCH THAT THE BOTTOM
      ! VELOCITY IS ZERO 

      allocate(eta_lcl(nobclsf)); eta_lcl=0.0_sp

      rhointcur= 0.0_sp
      rhointnxt = 0.0_sp

      DO i=nobclsf,1,-1  ! COUNT BACKWARDS
         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
         
!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
      

      DEALLOCATE(eta_lcl)

   END IF


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

   DO i=1,n

!
!--2 SOLID BOUNDARY EDGES------------------------------------------------------|
!
!Q&C< commented on 06/22/04
!   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

!Q&C< commented on 06/22/04
!           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

!           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
!           UAS(I)=UAF(I)
!           VAS(I)=VAF(I)

21        CONTINUE
          END IF
       END DO


!==============================================================================|
   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.

       do i= 1, n
         do k =1, kbm1
!Q&C< commented on 06/22/04
!          if(isbce(i).eq.3) then
!           uf(i,k)=0.0_SP
!           vf(i,k)=0.0_SP
!          end if
!Q&C> 06/22/04

          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

!Q&C< commented on 06/22/04
!          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

!          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
!          us(i,k)=uf(i,k)
!           vs(i,k)=vf(i,k)
!Q&C> 06/22/04

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
!Q&C< commented on 06/22/04
!          if(isbce(i).eq.3) then
!          u(i,k)=0.0_SP
!          v(i,k)=0.0_SP
!          end if
!Q&C> 06/22/04
          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

!Q&C< commented on 06/22/04
!          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

!          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
!          us(i,k)=u(i,k)
!          vs(i,k)=v(i,k)
!Q&C> 06/22/04

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
        
        
        wusurf = -wusurf * ramp * 0.001_sp
        wvsurf = -wvsurf * ramp * 0.001_sp
        
     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)
     END IF




!
!-- Set Groundwater flux ------------------------------------------------------|
!
     
     IF (groundwater_on) THEN
        CALL update_groundwater(inttime,bfwdis,gw_temp=bfwtmp,gw_salt=bfwslt)
        bfwdis = ramp * bfwdis
     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)
   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_gcy"

Here is the call graph for this function: