BODYFORCE-transform acceleration to stress

[zhangtianhao]

Hi there,

When using BODYFORCE, we can apply the surface stress to the whole water column. And BODYFORCE requirs sustr and svstr which can be set in smflux.F.

However, as for me, I only want to set a varying horizontal acceration in the whole water column, so I don't have sustr/svstr data but I obtain acceleration ax and ay data.

Now I modify the smflux.F as followed:

Code: Select all

!!!!!!!Hoty modified at 0614!!!!!!

      open(55,file='/ubuntu/ROMS/Projects/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*(abs(h(i-1,j))+abs(h(i,j))+z_w(i,j,N(ng))+z_w(i-1,j,N(ng)))*ax(k)
	    END IF
	  END DO

In short, the formula is sustr (m2/s2) = acceleration (m/s2) * whole depth (m)

I wonder that if is right/ How do you deal with the acceleration when using BODYFORCE?

Thanks and look forward to your reply.

Hoty

[wilkin]

There is a different approach to imposing an effective horizontal acceleration that is uniform throughout the vertical water column. This is to introduce a horizontal gradient in the air pressure (e.g. within ana_pair.h) and impose this force via #define ATM_PRESS.

I use this method to impose a time varying oscillatory force in small domain doubly periodic (i.e. effectively 1-D vertical) configurations such as BIO_TOY.

The air pressure sets the surface boundary condition to the hydrostatic relation, but it's role in the momentum equation is only via its gradient ... on the RHS the term 1/rho0 grad(p).

Creating an artificial P_air field that is a linear function of, say, x, such as P_air = ax will introduce a uniform acceleration 1/rho0 dp/dx = a/rho0. The units of a are Pascals/m and you can verify that a/rho0 is m/s^2.

Be careful if you are using P_air for anything else, like BULK_FLUXES, but it sounds as though you are not.

If you attempt this, be aware that ana_pair.h is working with air pressure in units of millibars so you need to convert Pa to millibars by multiplying by 0.01.

If you do take this approach, and you are doing so in a periodic channel, you need to modify ana_pair.h NOT to apply the periodic boundary conditions exchange, which is the block of code at the end of the Functional:

Code: Select all

#if !defined USE_PAIR_TO_APPLY_ACCELERATION
/* Do not do the periodic boundary condition exchange  
   or it breaks the constant x-direction pressure forcing */
      IF (EWperiodic(ng).or.NSperiodic(ng)) THEN 
       CALL exchange_r2d_tile (ng, tile,                               &
     &                          LBi, UBi, LBj, UBj,                    &
     &                          Pair)
      END IF
#endif

My code for this is messy being a special case so I won't post it here, but if you want to see how I got this to work message me directly and I'll share my clunky ana_pair.h

[zhangtianhao]

There is a different approach to imposing an effective horizontal acceleration that is uniform throughout the vertical water column. This is to introduce a horizontal gradient in the air pressure (e.g. within ana_pair.h) and impose this force via #define ATM_PRESS.

I use this method to impose a time varying oscillatory force in small domain doubly periodic (i.e. effectively 1-D vertical) configurations such as BIO_TOY.

The air pressure sets the surface boundary condition to the hydrostatic relation, but it's role in the momentum equation is only via its gradient ... on the RHS the term 1/rho0 grad(p).

Creating an artificial P_air field that is a linear function of, say, x, such as P_air = ax will introduce a uniform acceleration 1/rho0 dp/dx = a/rho0. The units of a are Pascals/m and you can verify that a/rho0 is m/s^2.

Be careful if you are using P_air for anything else, like BULK_FLUXES, but it sounds as though you are not.

If you attempt this, be aware that ana_pair.h is working with air pressure in units of millibars so you need to convert Pa to millibars by multiplying by 0.01.

If you do take this approach, and you are doing so in a periodic channel, you need to modify ana_pair.h NOT to apply the periodic boundary conditions exchange, which is the block of code at the end of the Functional:

Code: Select all

#if !defined USE_PAIR_TO_APPLY_ACCELERATION
/* Do not do the periodic boundary condition exchange  
   or it breaks the constant x-direction pressure forcing */
      IF (EWperiodic(ng).or.NSperiodic(ng)) THEN 
       CALL exchange_r2d_tile (ng, tile,                               &
     &                          LBi, UBi, LBj, UBj,                    &
     &                          Pair)
      END IF
#endif

My code for this is messy being a special case so I won't post it here, but if you want to see how I got this to work message me directly and I'll share my clunky ana_pair.h

Thanks,

I will take a try and then post my results.