the defination of density rho

[zhangtianhao]

Hi there,

I have a problem and maybe it's easy to you: I want to modify the code and use the density of seawater in fourmula, should I use rho(:,:,:)+1000 or just rho? I wonder that if the rho is the density anomaly I don't find the defination of in forum.

BTW, I want to use the density to represent the mass but not use it in pressure terms.

Thanks

[kate]

The variable "rho" in the model is insitu density anomaly. You add rho0 to it to get full density. Note that rho0 may or may not be 1000.

[zhangtianhao]

The variable "rho" in the model is insitu density anomaly. You add rho0 to it to get full density. Note that rho0 may or may not be 1000.

Thanks, Kate, rho0 equals 1025 in my case, but I find in rho_eos.F file:

Code: Select all

den(i,k)=den1(i,k)*bulk(i,k)*cff
den(i,k)=den(i,k)-1000.0_r8

It seems that in this code file, they do not use den-rho0 but den-1000..

[kate]

OK, but I'm surprised.

[zhangtianhao]

Hi all,

I found that when I set a (rho+1000)*acceleration term in the equation, ROMS showed rho speed trouble, I wonder that how to check this problem.

Thanks

[arango]

Are you trying to change ROMS equation of state? The definition of in situ density anomaly or potential density anomaly is that a1000 kg m-3 value is subtracted to the density. The in situ density anomaly, as a field, is only used in ROMS in the pressure gradient algorithm where its value is not essential. We are only interested in its horizontal gradient! But, you cannot change the ROMS equation of state without looking at how the pressure gradient algorithm is formulated. Also, the potential density is used in the isopycnic diffusion along neutral surfaces in the Redi tensor.

It is well known from the beggigning of oceanography that we cannot measure the seawater density directly with a sufficient level of accuracy to determine its horizontal gradient. That's the reason that we measure temperature and salinity as a function of depth and compute density with the equation of state for seawater.

If you need the density without the pressure terms (secant bulk modulus) use pden, which is the density anomaly at standard one-atmosphere pressure.

The parameter rho0 is a different density value (not an anomaly value) that it is used in the Boussinesq approximation of the governing equations. The density gradients are ignored except when multiplied by gravity. Rho0 is also used when activating the linear equation of state (NONLIN_EOS is off). Its value is set-up in ROMS standard input file roms.in. If you preffer a rho0 value of 1000 kg m-3, do so. It is a user tunable parameter!

[zhangtianhao]

Are you trying to change ROMS equation of state? The definition of in situ density anomaly or potential density anomaly is that a1000 kg m-3 value is subtracted to the density. The in situ density anomaly, as a field, is only used in ROMS in the pressure gradient algorithm where its value is not essential. We are only interested in its horizontal gradient! But, you cannot change the ROMS equation of state without looking at how the pressure gradient algorithm is formulated. Also, the potential density is used in the isopycnic diffusion along neutral surfaces in the Redi tensor.

It is well known from the beggigning of oceanography that we cannot measure the seawater density directly with a sufficient level of accuracy to determine its horizontal gradient. That's the reason that we measure temperature and salinity as a function of depth and compute density with the equation of state for seawater.

If you need the density without the pressure terms (secant bulk modulus) use pden, which is the density anomaly at standard one-atmosphere pressure.

The parameter rho0 is a different density value (not an anomaly value) that it is used in the Boussinesq approximation of the governing equations. The density gradients are ignored except when multiplied by gravity. Rho0 is also used when activating the linear equation of state (NONLIN_EOS is off). Its value is set-up in ROMS standard input file roms.in. If you preffer a rho0 value of 1000 kg m-3, do so. It is a user tunable parameter!

Thank you very much for your detailed answer,

I want to modify momentum equation in rhs3d by adding an forcing term (a constant horizental acceleration a), and density is required in this term. According to this, should I use (pden+1000) in this term, or should I use rho0 because the acceleration is a body force, like g (Boussinesq approximation)?

Thanks again

Hoty

[wilkin]

Where ROMS expresses anything in kinematic units (per unit mass), such as converting stress to friction velocity (tau/rho0), or heat flux (W/m2) to degrees meters/s (Hf/(Cp*rho0) in the surface temperature boundary condition, it always uses the constant parameter rho0. Likewise, to calculate Brunt-Väisälä frequency the factor g/rho0 is used.

If you are trying to add a "horizontal acceleration", are you aware of the existing BODYFORCE option?

[zhangtianhao]

Where ROMS expresses anything in kinematic units (per unit mass), such as converting stress to friction velocity (tau/rho0), or heat flux (W/m2) to degrees meters/s (Hf/(Cp*rho0) in the surface temperature boundary condition, it always uses the constant parameter rho0. Likewise, to calculate Brunt-Väisälä frequency the factor g/rho0 is used.

If you are trying to add a "horizontal acceleration", are you aware of the existing BODYFORCE option?

Thanks for your reply, I checked BODYFORCE by your guidence, I found the code in rhs3d.F about BODYFORCE:

Code: Select all

# ifdef BODYFORCE
      DO j=JstrV-1,Jend
        DO i=IstrU-1,Iend
          wrk(i,j)=0.0_r8
        END DO
      END DO
      DO k=N(ng),levsfrc(ng),-1
        DO j=JstrV-1,Jend
          DO i=IstrU-1,Iend
            wrk(i,j)=wrk(i,j)+Hz(i,j,k)
          END DO
        END DO
      END DO
      DO j=Jstr,Jend
        DO i=IstrU,Iend
          cff=0.25_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
          cff1=1.0_r8/(cff*(wrk(i-1,j)+wrk(i,j)))
          Uwrk(i,j)=sustr(i,j)*cff1
        END DO
      END DO
      DO k=levsfrc(ng),N(ng)
        DO j=Jstr,Jend
          DO i=IstrU,Iend
            cff=Uwrk(i,j)*(Hz(i  ,j,k)+Hz(i-1,j,k))
            ru(i,j,k,nrhs)=ru(i,j,k,nrhs)+cff

It seems that the forcing term ru on the right hand side of the equation euqals τ*A*Hz/Wrk，which has a force (N) dimension (Pa*m2*m/m).

As for me, BODYFORCE requirs surface stress and transforms it to body force, BUT I input acceleration a(m/s2) (because I just have horizontal acceleration data, I let sustr=au and svstr=ay mandatorily ), not stress in smflux.F, and I modify the code in order to satisfy the dimensional homogeneity: NOTICE the second line from the bottom, ru= rho0*acceleration*A*Wrk*Hz/Wrk

Code: Select all

# ifdef BODYFORCE
      DO j=JstrV-1,Jend
        DO i=IstrU-1,Iend
          wrk(i,j)=0.0_r8
        END DO
      END DO
      DO k=N(ng),levsfrc(ng),-1
        DO j=JstrV-1,Jend
          DO i=IstrU-1,Iend
            wrk(i,j)=wrk(i,j)+Hz(i,j,k)
          END DO
        END DO
      END DO
      DO j=Jstr,Jend
        DO i=IstrU,Iend
          cff=0.25_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
          cff1=1.0_r8/(cff*(wrk(i-1,j)+wrk(i,j)))
          Uwrk(i,j)=sustr(i,j)*cff1[/b]
        END DO
      END DO
      DO k=levsfrc(ng),N(ng)
        DO j=Jstr,Jend
          DO i=IstrU,Iend
            cff=Uwrk(i,j)*(Hz(i  ,j,k)+Hz(i-1,j,k))*rho0*(wrk(i-1,j)+wrk(i,j))*0.5_r8   !!!!!!!!!!!!NOTICE THIS MODIFICATION
            ru(i,j,k,nrhs)=ru(i,j,k,nrhs)+cff

But, ROMS showed rho speed trouble. I wonder that if there is anything I misunderstood.

Thanks

Hoty

[wilkin]

It seems that the forcing term ru on the right hand side of the equation euqals τ*A*Hz/Wrk，which has a force (N) dimension (Pa*m2*m/m).

No, at this point in the code sustr is a kinematic stress tao/rho0 and has units m2/s2. It is not Pascals.

See in bulk_flux.F and ana_smflux.h that Tau gets divided by rho0 to carry kinematic stress internally in ROMS.

[zhangtianhao]

Thanks for all of your advice,

I use BODYFORCE now, and I correct the τ unit. But ROMS still blowed up.

In smflx.h, I set sustr=a*h, (a means acceleration, h means depth, unit=m/s2*m), so the sustr has kinetic unit, as shown below:

Code: Select all

      open(55,file='/ubuntu/COAWST/Projects/acceleration/EW01.txt',status='old')
      read(55,*)(ax(i),i=1,40001)    
      close(55)

      DO j=JstrT,JendT
        DO i=IstrP,IendT
	  DO k=1,40001
	    kk = REAL(k)
	    IF(iic(ng).eq.kk) THEN   
              sustr(i,j)=0.5_r8*(h(i-1,j)+h(i,j))*ax(k)
	    END IF
	  END DO

Thus I don't need to modify BODYFORCE sustr code, (But I modify sbstr code, I don't want to transform it to bodyforce, it is still used in 2d equation), as shown below:

Code: Select all

!  Apply bottom stress as a bodyforce: determine the thickness (m)
!  of the bottom layer; then add in bottom stress as a bodyfoce.
!
      DO j=JstrV-1,Jend
        DO i=IstrU-1,Iend
          wrk(i,j)=0.0_r8
        END DO
      END DO
      DO k=1,levbfrc(ng)
        DO j=JstrV-1,Jend
          DO i=IstrU-1,Iend
            wrk(i,j)=wrk(i,j)+Hz(i,j,k)
          END DO
        END DO
      END DO
      DO j=Jstr,Jend
        DO i=IstrU,Iend
          cff=0.25_r8*(pm (i-1,j)+pm (i,j))*                            &
     &                (pn (i-1,j)+pn (i,j))
          cff1=1.0_r8/(cff*(wrk(i-1,j)+wrk(i,j)))
          Uwrk(i,j)=bustr(i,j)*cff1
        END DO
      END DO
      DO k=1,levbfrc(ng)
        DO j=Jstr,Jend
          DO i=IstrU,Iend
            cff=Uwrk(i,j)*(Hz(i  ,j,k)+                                 &
     &                     Hz(i-1,j,k))
!  hoty modified at 2020-06-11!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            ru(i,j,k,nrhs)=ru(i,j,k,nrhs)

!!!!!!!!!NOTICE THAT I DONT ADD CFF TO RHS

But DIAG speed problem happened again. I wonder that if ther are something I still misunderstood. Maybe there are something wrong when I transform acceleration to τ.

[wilkin]

(But I modify sbstr code, I don't want to transform it to bodyforce, it is still used in 2d equation), as shown below:

All you need to do to prevent the bottom stress also being applied as a bodyforce is make LEVBFRC = 0 in roms.in:

Code: Select all

!------------------------------------------------------------------------------
! Body-force parameters. Used when CPP option BODYFORCE is activated.
!------------------------------------------------------------------------------
!
! LEVBFRC      Shallowest level to apply bottom momentum stress as a body-force.

[zhangtianhao]

(But I modify sbstr code, I don't want to transform it to bodyforce, it is still used in 2d equation), as shown below:

All you need to do to prevent the bottom stress also being applied as a bodyforce is make LEVBFRC = 0 in roms.in:

Code: Select all

!------------------------------------------------------------------------------
! Body-force parameters. Used when CPP option BODYFORCE is activated.
!------------------------------------------------------------------------------
!
! LEVBFRC      Shallowest level to apply bottom momentum stress as a body-force.

Thanks for your suggestions, it is really helpful!

But I still have two problems.

Firstly, it seems that when set LEVBFRC = 0 in roms.in, the bottom stress is still not activated in 2d momenten equation, in rhs3d.F( cff2= 0.0_r8):

Code: Select all

# ifdef BODYFORCE
          DO i=Istr,Iend
            cff=om_v(i,j)*on_v(i,j)
            cff1= 0.0_r8
            cff2= 0.0_r8
            rufrc(i,j)=rufrc(i,j)+cff1+cff2

Secondly, I still don't figure out the speed trouble problems. I think maybe there are some problems with my input stress in smflux.h (sustr=a*h), or this term:0.5_r8*(h(i-1,j)+h(i,j)) is not right? The code is shown below. I really wonder that how do you transform acceleration to sustr.

Code: Select all

      open(55,file='ubuntu/COAWST/Projects/EW/EW01.txt',status='old')
      read(55,*)(ax(i),i=1,40001)    
      close(55)

      DO j=JstrT,JendT
        DO i=IstrP,IendT
	  DO k=1,40001
	    kk = REAL(k)
	    IF(iic(ng).eq.kk) THEN   
              sustr(i,j)=0.5_r8*(h(i-1,j)+h(i,j))*ax(k)
	    END IF
	  END DO

Thanks again,
Hoty

[wilkin]

Code: Select all

    cff1= 0.0_r8
    cff2= 0.0_r8
    rufrc(i,j)=rufrc(i,j)+cff1+cff2

I can't find that code anywhere in my version of rhs3d.F.

[zhangtianhao]

Code: Select all

    cff1= 0.0_r8
    cff2= 0.0_r8
    rufrc(i,j)=rufrc(i,j)+cff1+cff2

I can't find that code anywhere in my version of rhs3d.F.

Sorry, I don't make myself clear,

Firstly , it seems that the miximun value of levbfrc is 1 but not 0,

Secondly, Line 2170 in rhs3d.F:

Code: Select all

# ifndef BODYFORCE
          DO i=Istr,Iend
            cff=om_v(i,j)*on_v(i,j)
            cff1= svstr(i,j)*cff
            cff2=-bvstr(i,j)*cff
            rvfrc(i,j)=rvfrc(i,j)+cff1+cff2

It seems that when undefine BODYFORCE cff2(bottem stress) is introduced into right hand side of 2d equation.
That means if I define BODYFORCE, bottem stress will not be added in rhs of 2d eq, but added in the lowest layer in 3d euqation.

I wonder that if these two method are equivalent (add in 2d eq or the lowest layer in 3d eq)

[wilkin]

I see your point. Then you might have to code some workaround, like ...

Code: Select all

# if !defined BODYFORCE
! the original code ...
          DO i=Istr,Iend
            cff=om_v(i,j)*on_v(i,j)
            cff1= svstr(i,j)*cff
            cff2=-bvstr(i,j)*cff
            rvfrc(i,j)=rvfrc(i,j)+cff1+cff2
      ...      
#else 
! the code you want, whatever that is
...

But before you do that think it through carefully. Applying the bottom stress as a bodyforce in the bottom layer (LEVBFRC=1) will be extremely close to, if not identical to, having a stress on the z=-h face of the lowest layer. Isn't that what you want?

[zhangtianhao]

I see your point. Then you might have to code some workaround, like ...

Code: Select all

# if !defined BODYFORCE
! the original code ...
          DO i=Istr,Iend
            cff=om_v(i,j)*on_v(i,j)
            cff1= svstr(i,j)*cff
            cff2=-bvstr(i,j)*cff
            rvfrc(i,j)=rvfrc(i,j)+cff1+cff2
      ...      
#else 
! the code you want, whatever that is
...

But before you do that think it through carefully. Applying the bottom stress as a bodyforce in the bottom layer (LEVBFRC=1) will be extremely close to, if not identical to, having a stress on the z=-h face of the lowest layer. Isn't that what you want?

Thanks, I will take a try and have a look at the difference between the two methods.

Thanks agian

Hoty