--- MITgcm/verification/hs94.cs-32x32x5/code/external_forcing.F	2001/06/06 20:31:48	1.4
+++ MITgcm/verification/hs94.cs-32x32x5/code/external_forcing.F	2001/07/06 22:13:37	1.5
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/verification/hs94.cs-32x32x5/code/Attic/external_forcing.F,v 1.4 2001/06/06 20:31:48 adcroft Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/verification/hs94.cs-32x32x5/code/Attic/external_forcing.F,v 1.5 2001/07/06 22:13:37 jmc Exp $
 C $Name:  $
 
 #include "CPP_OPTIONS.h"
@@ -46,18 +46,21 @@
 C     _RL dist2
 C     _RL decayFac
 C     _RL velDragHeightFac
-      _RL termP,kV,kF
+      _RL recip_P0g,termP,kV,kF,sigma_b
 
 C--   Forcing term(s)
       kF=1. _d 0/86400. _d 0
-      DO J=jMin,jMax
-       DO I=iMin,iMax
-        IF ( HFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN
-C        termP=0.5*( rF(kLev) + min( rF(kLev+1) ,
-C    &           min(H(I,J,bi,bj),H(I,J-1,bi,bj))            ) )
-         termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )
-         kV=kF*MAX(0. _d 0,
-     &   (termP*recip_Rcol(I,J,bi,bj)-0.7 _d 0)/(1. _d 0-0.7 _d 0) )
+      sigma_b = 0.7 _d 0
+c     DO J=jMin,jMax
+c      DO I=iMin,iMax
+      DO J=1,sNy
+       DO I=1,sNx+1
+        IF ( hFacW(i,j,kLev,bi,bj) .GT. 0. ) THEN
+         recip_P0g=MAX(recip_Rcol(I,J,bi,bj),recip_Rcol(I-1,J,bi,bj))
+         termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
+     &                   +rF(kLev+1)*recip_P0g )
+c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g
+         kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
          gU(i,j,kLev,bi,bj)=gU(i,j,kLev,bi,bj)
      &                      -kV*uVel(i,j,kLev,bi,bj)
         ENDIF
@@ -109,18 +112,21 @@
 C     _RL dist2
 C     _RL decayFac
 C     _RL velDragHeightFac
-      _RL termP,kV,kF
+      _RL recip_P0g,termP,kV,kF,sigma_b
 
 C--   Forcing term(s)
       kF=1. _d 0/86400. _d 0
-      DO J=jMin,jMax
-       DO I=iMin,iMax
-        IF ( HFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN
-C        termP=0.5*( rF(kLev) + min( rF(kLev+1) ,
-C    &           min(H(I,J,bi,bj),H(I,J-1,bi,bj))            ) )
-         termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )
-         kV=kF*MAX(0. _d 0,
-     &   (termP*recip_Rcol(I,J,bi,bj)-0.7 _d 0)/(1. _d 0-0.7 _d 0) )
+      sigma_b = 0.7 _d 0
+c     DO J=jMin,jMax
+c      DO I=iMin,iMax
+      DO J=1,sNy+1
+       DO I=1,sNx
+        IF ( hFacS(i,j,kLev,bi,bj) .GT. 0. ) THEN
+         recip_P0g=MAX(recip_Rcol(I,J,bi,bj),recip_Rcol(I,J-1,bi,bj))
+         termP=0.5 _d 0*( MIN(rF(kLev)*recip_P0g,1. _d 0)
+     &                   +rF(kLev+1)*recip_P0g )
+c        termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )*recip_P0g
+         kV=kF*MAX( 0. _d 0, (termP-sigma_b)/(1. _d 0-sigma_b) )
          gV(i,j,kLev,bi,bj)=gV(i,j,kLev,bi,bj)
      &                      -kV*vVel(i,j,kLev,bi,bj)
         ENDIF
@@ -164,26 +170,26 @@
 C     == Local variables ==
 C     Loop counters
       INTEGER I, J
-      _RL thetaLim,kT,ka,ks,term1,term2,thetaEq,termP,rSurf
+      _RL thetaLim,kT,ka,ks,sigma_b,term1,term2,thetaEq,termP
 
 C--   Forcing term(s)
-      rSurf=1. _d 5
       ka=1. _d 0/(40. _d 0*86400. _d 0)
       ks=1. _d 0/(4. _d 0 *86400. _d 0)
+      sigma_b = 0.7 _d 0
       DO J=jMin,jMax
        DO I=iMin,iMax
-       term1=60. _d 0*(sin(yC(I,J,bi,bj)*deg2rad)**2)
-C      termP=0.5*( rF(kLev) + min( rF(kLev+1) , H(I,J,bi,bj) ) )
-       termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )
-       term2=10. _d 0*log(termP/rSurf)
-     &          *(cos(yC(I,J,bi,bj)*deg2rad)**2)
-       thetaLim = 200. _d 0/ ((termP/rSurf)**(2. _d 0/7. _d 0))
-       thetaEq=315. _d 0-term1-term2
-       thetaEq=MAX(thetaLim,thetaEq)
-        kT=ka+(ks-ka)
-     &    *MAX(0. _d 0,
-     &   (termP*recip_Rcol(I,J,bi,bj)-0.7 _d 0)/(1. _d 0-0.7 _d 0) )
-     &    *COS((yC(I,J,bi,bj)*deg2rad))**4
+         term1=60. _d 0*(sin(yC(I,J,bi,bj)*deg2rad)**2)
+         termP=0.5 _d 0*( rF(kLev) + rF(kLev+1) )
+         term2=10. _d 0*log(termP/atm_po)
+     &            *(cos(yC(I,J,bi,bj)*deg2rad)**2)
+         thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
+         thetaEq=315. _d 0-term1-term2
+         thetaEq=MAX(thetaLim,thetaEq)
+         termP=0.5 _d 0*( MIN(rF(kLev),Ro_surf(I,J,bi,bj))+rF(kLev+1) ) 
+         kT=ka+(ks-ka)
+     &     *MAX(0. _d 0,
+     &       (termP*recip_Rcol(I,J,bi,bj)-sigma_b)/(1. _d 0-sigma_b) )
+     &     *COS((yC(I,J,bi,bj)*deg2rad))**4
          gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
      &        - kT*( theta(I,J,kLev,bi,bj)-thetaEq )
      &            *maskC(i,j,kLev,bi,bj)