C $Header: /u/gcmpack/MITgcm/verification/dome/code/obcs_calc.F,v 1.9 2011/05/24 20:31:33 jmc Exp $ C $Name: $ #include "OBCS_OPTIONS.h" SUBROUTINE OBCS_CALC( futureTime, futureIter, & uVel, vVel, wVel, theta, salt, & myThid ) C *==========================================================* C | SUBROUTINE OBCS_CALC C | o Calculate future boundary data at open boundaries C | at time = futureTime C *==========================================================* IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "OBCS_PARAMS.h" #include "OBCS_GRID.h" #include "OBCS_FIELDS.h" #include "EOS.h" C == Routine arguments == INTEGER futureIter _RL futureTime _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) INTEGER myThid #ifdef ALLOW_OBCS C == Local variables == INTEGER bi, bj INTEGER I, J ,K _RL rampTime2 _RL rexpt _RL hinit, delh _RL z(nr) _RL Dmax,Dinf,dTemp,gp_inflow,Lrho _RL Width,x,Xcenter,Fz,zt,Rit C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC',myThid) #endif c Total depth Dmax = 600. _d 0 c width of flow transition delh=5. _d 0 c total depth of active layer Dinf=300. _d 0 c g_prime of inflow dTemp=(2. _d 0/rhonil)/(2. _d -4) gp_inflow=tAlpha*gravity*dTemp c Deformation radius Lrho=sqrt(gp_inflow*Dinf)/f0 c dimensional width Width=100. _d 3 c nondimensional width Width=Width/Lrho c print *,'Width=',Width c coordinate of center of embayment Xcenter = 1700. _d 3 c Critical Richardson number Rit = 1. _d 0/3. _d 0 c Ttarget=-5.0e6 c areaCh=Width*Lrho*Dmax c relaxtime=180000.0 c relaxConst=1.0/(areaCh*relaxtime) z(1) = -delR(1)/2. _d 0 do j=2,nr Caja NOTE: This gives the depths of W points, not U,V,T,S,p points C USe instead drC but requires different indexing. See documentation C in AJAs head z(j) = z(j-1) - delR(j) enddo c ramptime for velocity rampTime2 = 4. _d 0*44567. _d 0 rexpt=1. _d 0/exp(futureTime/rampTime2) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C Eastern OB IF (useOrlanskiEast) THEN CALL ORLANSKI_EAST( & bi, bj, futureTime, & uVel, vVel, wVel, theta, salt, & myThid ) ELSE DO K=1,Nr DO J=1-Oly,sNy+Oly OBEu(J,K,bi,bj)=0. c & +Uinflow c & *0.5*(1.0 - (exp((z(k)-hinit)/delh) c & - exp(-(z(k)-hinit)/delh)) c & /(exp((z(k)-hinit)/delh) c & + exp(-(z(k)-hinit)/delh))) c &*(1.0 - rexpt*rexpt)/(1.0 + rexpt*rexpt) OBEv(J,K,bi,bj)=0. OBEt(J,K,bi,bj)=tRef(K) c & +Tinflow c & *0.5*(1.0 - (exp((z(k)-hinit)/delh) c & - exp(-(z(k)-hinit)/delh)) c & /(exp((z(k)-hinit)/delh) c & + exp(-(z(k)-hinit)/delh))) c &*(1.0 - rexpt*rexpt)/(1.0 + rexpt*rexpt) OBEs(J,K,bi,bj)=sRef(K) #ifdef ALLOW_NONHYDROSTATIC OBEw(J,K,bi,bj)=0. #endif ENDDO ENDDO ENDIF C Western OB IF (useOrlanskiWest) THEN CALL ORLANSKI_WEST( & bi, bj, futureTime, & uVel, vVel, wVel, theta, salt, & myThid ) ELSE DO K=1,Nr DO J=1-Oly,sNy+Oly OBWu(J,K,bi,bj)=0. OBWv(J,K,bi,bj)=0. OBWt(J,K,bi,bj)= Tref(k) #ifdef ALLOW_NONHYDROSTATIC OBWw(J,K,bi,bj)=0.0 #endif ENDDO ENDDO ENDIF C Northern OB, template for forcing IF (useOrlanskiNorth) THEN CALL ORLANSKI_NORTH( & bi, bj, futureTime, & uVel, vVel, wVel, theta, salt, & myThid ) ELSE DO K=1,Nr DO I=1-Olx,sNx+Olx c Make center of embayment x=0 x=(xC(I,1,bi,bj)-Xcenter)/Lrho +Width/2. _d 0 IF ((x.GE.0.).AND.(x.LE.Width)) THEN hinit=Dinf*(exp(-x)) - Dmax c if (k.eq.5) then c print *,'x=',x,' xC',xC(I,1,bi,bj),' I=',I c print *,'hinit=',hinit c endif zt=(z(k) + Dmax - (hinit + Dmax))/(hinit+Dmax) IF (zt.GE.(Rit/(2. _d 0 -Rit))) THEN Fz=1. ELSE IF (zt.GE.(-Rit/(2.+Rit))) THEN Fz=1. _d 0/Rit*zt/(zt+1. _d 0) + 0.5 _d 0 ELSE Fz=0. ENDIF ENDIF ELSE Fz=1. ENDIF c if ((x.ge.0).and.(x.le.Width)) then c print *,'z(k)',z(k),'zt',zt,'Fz',Fz c endif OBNv(I,K,bi,bj)=0. & - sqrt(gp_inflow*Dinf)*exp(-x) & *(1. _d 0 - Fz) c if ((x.ge.0).and.(x.le.Width)) then c if (k.eq.5) then c print *,'V=',OBNv(I,K,bi,bj) c endif c endif OBNu(I,K,bi,bj)=0. IF (tRef(K).LE. (-dTemp*(1. _d 0 - Fz))) THEN OBNt(I,K,bi,bj) = tRef(K) ELSE OBNt(I,K,bi,bj) = -dTemp*(1. _d 0 - Fz) ENDIF c if ((x.ge.0).and.(x.le.Width)) then c print *,'T=',OBNt(I,K,bi,bj) c endif c OBNt(I,K,bi,bj)=tRef(K) c & - dTemp c & *(1.0 - Fz) OBNs(I,K,bi,bj)=sRef(K) + 1. _d 0*(1. _d 0 - Fz) #ifdef ALLOW_NONHYDROSTATIC OBNw(I,K,bi,bj)=0. #endif ENDDO ENDDO ENDIF C Southern OB, template for forcing IF (useOrlanskiSouth) THEN CALL ORLANSKI_SOUTH( & bi, bj, futureTime, & uVel, vVel, wVel, theta, salt, & myThid ) ELSE DO K=1,Nr DO I=1-Olx,sNx+Olx OBSu(I,K,bi,bj)=0. OBSv(I,K,bi,bj)=0. OBSt(I,K,bi,bj)=tRef(K) OBSs(I,K,bi,bj)=sRef(K) #ifdef ALLOW_NONHYDROSTATIC OBSw(I,K,bi,bj)=0. #endif ENDDO ENDDO ENDIF C-- end bi,bj loops. ENDDO ENDDO #ifdef ALLOW_OBCS_BALANCE IF ( useOBCSbalance ) THEN CALL OBCS_BALANCE_FLOW( futureTime, futureIter, myThid ) ENDIF #endif /* ALLOW_OBCS_BALANCE */ #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC',myThid) #endif #endif /* ALLOW_OBCS */ RETURN END