[MITgcm-support] the effective number of bits in obcs_calc.F

[H.W. Chou]

Dear all,

Good evening, this is Shawn from Hokkaido University. 
For setting the K1 tide in high latitude area(53 degree), 
I followed the equation of tidal residual flow :
   u=-(c/H)*(eta-eta0*sin(omega*time)) in east open boundary
and 
   u=(c/H)*eta  as radiation condition in west open boundary,
which omega=2*pi/(period of tide) rad s-1  is frequency of tide
c=((gH(omega)^2)/((omega)^2-f^2))^(1/2), is the linear long wave speed
H is the deep water depth,
eta is the sea-surface elevation 
eta0 = u0*H/c where I set u0 is 5.0 cm s-1

I add a part of code in program obcs_calc.F as below:

In eastern part: 

#ifdef ALLOW_OBCS_EAST
C     Eastern OB
#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: East',myThid)
#endif
      IF (useOrlanskiEast) THEN
#ifdef ALLOW_ORLANSKI
        CALL ORLANSKI_EAST(
     &          bi, bj, futureTime,
     &          uVel, vVel, wVel, theta, salt,
     &          myThid )
#endif
#ifdef ALLOW_NEST_CHILD
      ELSEIF ( useNEST_CHILD ) THEN
        DO k=1,Nr
          DO j=1-OLy,sNy+OLy
            IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
              OBEu(j,k,bi,bj)= U_F1(j,k,2)
              OBEv(j,k,bi,bj)= V_F1(j,k,2)
              OBEt(j,k,bi,bj)= T_F1(j,k,2)
              OBEs(j,k,bi,bj)= S_F1(j,k,2)
#ifdef NONLIN_FRSURF
              OBEeta(j,bi,bj)= ETA_F1(j,1,2)
#endif
            ENDIF
          ENDDO
        ENDDO
#endif /* ALLOW_NEST_CHILD */
      ELSE
        DO k=1,Nr
          DO j=1-OLy,sNy+OLy
            I_obc = OB_Ie(j,bi,bj)
            IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
              IF (futureTime>=0.0) THEN
              zf = fCori(I_obc,j,bi,bj)
              zc = sqrt(gravity*zH*zomega2/(zomega2-zf*zf))
              zeta0 = zramp * zampl * (zH/zc)
              zeta = etaN(I_obc,j,bi,bj)
              OBEu(j,k,bi,bj)=(zc/zH)*
     &                    (zeta-zeta0*sin(zomega*futureTime)) 
              OBEv(j,k,bi,bj)=0.
              OBEt(j,k,bi,bj)=tRef(k)
              OBEs(j,k,bi,bj)=sRef(k)
#ifdef ALLOW_NONHYDROSTATIC
              OBEw(j,k,bi,bj)=0.
#endif
#ifdef NONLIN_FRSURF
              OBEeta(j,bi,bj)=0.
#endif
             ELSE
              OBEu(j,k,bi,bj)=0.
              OBEv(j,k,bi,bj)=0.
              OBEt(j,k,bi,bj)=tRef(k)
              OBEs(j,k,bi,bj)=sRef(k)
#ifdef ALLOW_NONHYDROSTATIC
              OBEw(j,k,bi,bj)=0.
#endif
#ifdef NONLIN_FRSURF
              OBEeta(j,bi,bj)=0.
#endif
             ENDIF 
            ENDIF
          ENDDO
        ENDDO
      ENDIF
#endif /* ALLOW_OBCS_EAST */

and the western part:
 #ifdef ALLOW_OBCS_WEST
C     Western OB
#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: West',myThid)
#endif
      IF (useOrlanskiWest) THEN
#ifdef ALLOW_ORLANSKI
        CALL ORLANSKI_WEST(
     &          bi, bj, futureTime,
     &          uVel, vVel, wVel, theta, salt,
     &          myThid )
#endif
#ifdef ALLOW_NEST_CHILD
      ELSEIF ( useNEST_CHILD ) THEN
        DO k=1,Nr
          DO j=1-OLy,sNy+OLy
            IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
              OBWu(j,k,bi,bj)= U_F1(j,k,1)
              OBWv(j,k,bi,bj)= V_F1(j,k,1)
              OBWt(j,k,bi,bj)= T_F1(j,k,1)
              OBWs(j,k,bi,bj)= S_F1(j,k,1)
#ifdef NONLIN_FRSURF
              OBWeta(j,bi,bj)= ETA_F1(j,1,1)
#endif
           ENDIF
          ENDDO
        ENDDO
#endif /* ALLOW_NEST_CHILD */
      ELSE
        DO k=1,Nr
          DO j=1-OLy,sNy+OLy
             I_obc = OB_Iw(j,bi,bj)
            IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
              IF (futureTime>0.0) THEN
              zf = fCori(I_obc,j,bi,bj)
              zc = sqrt(gravity*zH*zomega2/(zomega2-zf*zf))
              zeta0 = zramp * zampl * (zH/zc)
              zeta = etaN(I_obc,j,bi,bj) 
              OBWu(j,k,bi,bj) = (zc/zH)*etaN(I_obc+1,j,bi,bj) 
              OBWv(j,k,bi,bj)=0.
              OBWt(j,k,bi,bj)=tRef(k)
              OBWs(j,k,bi,bj)=sRef(k)
#ifdef ALLOW_NONHYDROSTATIC
              OBWw(j,k,bi,bj)=0.
#endif
#ifdef NONLIN_FRSURF
              OBWeta(j,bi,bj)=0.
#endif
            ELSE
              OBWu(j,k,bi,bj)=0.
              OBWv(j,k,bi,bj)=0.
              OBWt(j,k,bi,bj)=tRef(k)
              OBWs(j,k,bi,bj)=sRef(k)
#ifdef ALLOW_NONHYDROSTATIC
              OBWw(j,k,bi,bj)=0.
#endif
#ifdef NONLIN_FRSURF
              OBWeta(j,bi,bj)=0.
#endif
            ENDIF 
           ENDIF
          ENDDO
        ENDDO
      ENDIF
#endif /* ALLOW_OBCS_WEST */

The problem is, 
if I set the tide as M2 tide as few variables below, there is no problem 
in running:
      zomega = 1.405257046694307d-04
      zampl = 50.0d-2
      zomega2 = zomega*zomega
      zH = abs(rF(Nr+1))
      zramp = 1. d 0
However, if I changed the zomega into K1 tides frequency, all the 
results would be NaN directly:
      zomega = 0.728738727346275d-04  
      zampl = 50.0d-2
      zomega2 = zomega*zomega
      zH = abs(rF(Nr+1))
      zramp = 1. d 0
Is that because the zomega number too small? 

Thank you for any tips or hints. 
Sincerely,

Shawn