[MITgcm-support] the effective number of bits in obcs_calc.F
Martin Losch
Martin.Losch at awi.de
Fri Aug 17 09:19:46 EDT 2018
Hi Shawn,
did you get a reply?
I don’t see anything wrong in your code and I would expect that it should work independently of your zomega as long as your time step is small enough to resolve the tidal cycles properly (but you are going from shorter to longer period, so that shouldn’t be your problem). The difference between your zomegas is only about factor 1/2, so there is no adhoc reason why it should not work for the K1 frequency.
To help I would need you to generate a (small!!!) example setup, where you can reproduce your problem. If you can post the code and input directories somewhere for me to download (please don’t send them to this list), I could try to reproduce the problem.
Martin
> On 9. Aug 2018, at 07:43, H.W. Chou <shawn_chou at ees.hokudai.ac.jp> wrote:
>
>
> 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
>
>
>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
More information about the MITgcm-support
mailing list