which index of zeta should I use to add POT_TIDES

[xiaozhu557]

I am trying to add the POT_TIDES into set_tides.f90.
In my opinion, after the codes, as follow, of calculating the Potential Tides,

Code: Select all

          DO j=JstrR,JendR
            DO i=IstrR,IendR
#  ifdef TIDES_ASTRO
!  Convert Vu_sat from cycles to radians
              Ptide(i,j)=Ptide(i,j)+                                    &
     &                   ramp*POT_Tamp(i,j,itide)*                      &
     &                   COS(-POT_Tphase(i,j,itide)+                    &
     &                   2.0_r8*pi*Vu_sat(i,j,itide))
#  else
              Ptide(i,j)=Ptide(i,j)+                                    &
     &                   ramp*POT_Tamp(i,j,itide)*                      &
     &                   COS(omega-POT_Tphase(i,j,itide))
#  endif
#  ifdef MASKING
              Ptide(i,j)=Ptide(i,j)*rmask(i,j)
#  endif
            END DO
          END DO

we should add the variable Ptide to OCEAN(ng)%zeta, like

Code: Select all

OCEAN(ng)%zeta = OCEAN(ng)%zeta + Ptide

[/b].

But for the index i of OCEAN(ng)%zeta(:,:,i), it has three options, krhs, kstp, knew. Which one should I use?

Thanks.

[kate]

Right now, Ptide is added to the pressure P computed in prsgrd32.h. If you want to add it to zeta, you need to stop adding it to P in prsgrd32.h.

Looking at the top of step2d*.h:

Code: Select all

Drhs(i,j)=zeta(i,j,krhs)+h(i,j)

This tells you that krhs is what you want to add it to. But wait, just before the call to step2d, krhs is modified:

Code: Select all

                  IF (FIRST_2D_STEP) THEN
                    kstp(ng)=indx1(ng)
                  ELSE
                    kstp(ng)=3-indx1(ng)
                  END IF
                  knew(ng)=3
                  krhs(ng)=indx1(ng)

So you want whichever one was set to indx1 for the call to step2d before the call to set_tides:

Code: Select all

                  knew(ng)=next_indx1
                  kstp(ng)=3-knew(ng)
                  krhs(ng)=3
                  IF (iif(ng).lt.(nfast(ng)+1)) indx1(ng)=next_indx1

I'm thinking it should be kstp, but it might be best to check in a debugger or with print statements with this sort of logic.

[xiaozhu557]

Hello, Kate, thanks for your nice suggestions.

(1) For the original version, add Ptide to the P.
I found that if I define POT_TIDES, I have to define DJ_GRADPS. Because you have not add it in other files, prsgrd31.h, prsgrd40/2/4.h. I want to know why you did not add it in those files? If I want, could I simply add it to those files follow as the file prsgrd32.h?

(2) For the new version, what I want to change, I am still a little confused about your words "But wait, just .....". Because the step2d had been called two times, one for predictor step, the other for corrector step. Which one should I follow for the call to set_tides? Could I just simply copy one of the codes shown in your last post before the call to set_tides? Which one should I use?
Or whether I just use krhs is ok? Because I did not found the krhs was modified before the call to other subroutines, such set_vbc, set_zeta. They uses krhs or 1/2.

So could you explain more details or tell me how should I write the codes directly?

Thanks.

[wilkin]

The tidal potential force must be added to the pressure force because it's just that - a force. Its representation as an equilibrium tidal elevation is a convenient representation of the scalar field whose gradient has the same force effect, but it is not a time varying additive component of sea level

If you simply add the equilibrium tide directly to the the dynamic zeta you violate the continuity equation. You are adding and removing mass all over the ocean without any corresponding adjustment to the divergence of the depth integrated transport.

[xiaozhu557]

Thanks,Wilkin. I understand it clearly. Could you continue to answer my frst question in my last post?

(1) For the original version, add Ptide to the P.
I found that if I define POT_TIDES, I have to define DJ_GRADPS. Because you have not add it in other files, prsgrd31.h, prsgrd40/2/4.h. I want to know why you did not add it in those files? If I want, could I simply add it to those files follow as the file prsgrd32.h?

[wilkin]

It is my understanding that the DJ_GRADPS option is the recommended algorithm described in Section 5 "PGF" in Shchepetkin and McWilliams (2009):

Shchepetkin, A. F., and J. C. McWilliams, 2009: Computational Kernel Algorithms for Fine-Scale, Multiprocess, Longtime Oceanic Simulations. Handbook of Numerical Analysis, Computational Methods for the Atmosphere and the Oceans, 14, 121-183.

This paper is an excellent review and synthesis of the various algorithmic choices made in constructing the ROMS computational kernel.

The DJ_GRADPS algorithm is "robustly tolerant of hydrostatically inconsistent grids" and has "smaller truncation error and, more importantly, slower error growth in time for the flat stratification test cases."

To the best of my knowledge there are no ROMS users working on realistic applications that choose any option other than DJ_GRADPS.

[xiaozhu557]

Thank you so much Wilkin.