apply the atmosphere pressure in Storm simulation

[c4i]

Hi colleagues,

I am trying to run the ROMS model for a storm surge simulation in the Irish Sea.
How can I use the atmospheric pressure to force the free surface?
Has anybody done this before?
Any suggestion are welcome.

Thank you.

[kate]

I have been merging changes from Paul Budgell's code and noticed that he is doing this. He changed prsgrd32.h to include:

#ifdef SLP_GRAD
P(i,j,N(ng)) = P(i,j,N(ng)) + 100._r8*Pair(i,j)/rho0
#endif

This code goes after P(i,j,N(ng)) is set, before it is used to integrate downwards. Note that you need to pass Pair in through the arguments list along with the rest and that Pair is declared in mod_forces. Note also that I haven't tested this code.

[nobuhitomori]

Dear Kate,

I tried to include atmospheric pressure force to the free surface for storm surge in nearshore region.

I put the atmospheric pressure contribution just after the following sentence

Code: Select all

          P(i,j,N(ng))=GRho0*z_w(i,j,N(ng))+                            &
     &                 GRho*(rho(i,j,N(ng))+cff2)*                      &
     &                 (z_w(i,j,N(ng))-z_r(i,j,N(ng)))

as

Code: Select all

 P(i,j,N(ng)) = P(i,j,N(ng)) + 100._r8*Pair(i,j)/rho0

at "preliminary step" in prsgrd32.h of DJ_GRADP following your previous posted comments. I also changed the source codes of mod_forces and etc. I'm pretty sure that the Pair(i,j) is decreased from 1013hPa to 950hPa in this part but the free surface doesn't change in my test computation.

Do you have any comments on it?

I also have another question about pressure correction term of
P(i,j,N(ng)) = P(i,j,N(ng)) + 100._r8*Pair(i,j)/rho0
What's meaning of 100._r8?

[kate]

As I said, I haven't tested it, but I don't know why it wouldn't work. Can you check things out in a debugger?

For the other, ROMS wants to read Pair in millibars. This looks like it would convert Pair to Pascals for this calculation.

[azhang]

For my understanding, you need activate CPP Option "ATM_PRESS" so that ROMS declares arrays of PairG and Pair and also read in air pressure from forcing file if ANA_PAIR is not activated.

Since Pair and PairG are declared in Module mod_forces, I think the statement of "USE mod_forces" should be added in subroutine prsgrd*.h so the arrays Pair and PariG will be passed.

Good Luck

[nobuhitomori]

Dear kate and zhang

Thank you for valuable information. I didn't know the unit
difference between Pair and P. The ATM_PRESS is useful because
it looks better than using BULK_FLUXES.

I will play this modification and check.

[rongzr]

According to my experience, you need to edited this two files
(1)mode_forces.F--------Define Pair and PairG so that ROMS can read it
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!edited by Zengrui Rong 06/02/2008
#ifndef SOLVE3D
# if defined ATM_PRESS
real(r8), pointer :: Pair(:,:)
# ifndef ANA_PAIR
real(r8), pointer :: PairG(:,:,:)
# endif
# endif
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
(2) step2d_LF_AM3.h------- add the Air pressure term to the pressure gradient term
!=======================================================================
! Compute right-hand-side for the 2D momentum equations.
!=======================================================================
!
!-----------------------------------------------------------------------
! Compute pressure gradient terms.
!-----------------------------------------------------------------------
!
cff1=0.5_r8*g
cff2=1.0_r8/3.0_r8
DO j=Jstr,Jend
DO i=IstrU,Iend
rhs_ubar(i,j)=cff1*on_u(i,j)* &
& ((h(i-1,j)+ &
& h(i ,j))* &
& (gzeta(i-1,j)- &
& gzeta(i ,j))+ &
# if defined VAR_RHO_2D && defined SOLVE3D
& (h(i-1,j)- &
& h(i ,j))* &
& (gzetaSA(i-1,j)+ &
& gzetaSA(i ,j)+ &
& cff2*(rhoA(i-1,j)- &
& rhoA(i ,j))* &
& (zwrk(i-1,j)- &
& zwrk(i ,j)))+ &
# endif
& (gzeta2(i-1,j)- &
& gzeta2(i ,j)))+ &
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! impose atmospheric pressure onto sea surface
!!!!! added by Zengrui Rong 05/30/2008
# ifndef SOLVE3D
# ifdef ATM_PRESS
& 0.5_r8*on_u(i,j)*100.0/rho0* &
& ( h(i-1,j)+h(i,j)+ &
& gzeta(i-1,j)+gzeta(i,j))* &
& (Pair(i-1,j)-Pair(i,j))+ &
# endif
# endif
& 0.0_r8
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# ifdef DIAGNOSTICS_UV
DiaU2rhs(i,j,M2pgrd)=rhs_ubar(i,j)
# endif
END DO
IF (j.ge.JstrV) THEN
DO i=Istr,Iend
rhs_vbar(i,j)=cff1*om_v(i,j)* &
& ((h(i,j-1)+ &
& h(i,j ))* &
& (gzeta(i,j-1)- &
& gzeta(i,j ))+ &
# if defined VAR_RHO_2D && defined SOLVE3D
& (h(i,j-1)- &
& h(i,j ))* &
& (gzetaSA(i,j-1)+ &
& gzetaSA(i,j )+ &
& cff2*(rhoA(i,j-1)- &
& rhoA(i,j ))* &
& (zwrk(i,j-1)- &
& zwrk(i,j )))+ &
# endif
& (gzeta2(i,j-1)- &
& gzeta2(i,j )))+ &
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! impose atmospheric pressure onto sea surface
!!!!! added by Zengrui Rong 05/30/2008
# ifndef SOLVE3D
# ifdef ATM_PRESS
& 0.5_r8*om_v(i,j)*100.0/rho0* &
& ( h(i,j-1)+h(i,j)+ &
& gzeta(i,j-1)+gzeta(i,j))* &
& (Pair(i,j-1)-Pair(i,j)) + &
# endif
# endif
& 0.0_r8
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!




Note that this only works for the 2D case.

[kate]

It's true that Paul's method only works in 3-D. Also, I have bulk_flux on because that's how I compute my stresses, so I have Pair available to me anyway. We don't seem to be able to foresee all the options that people need, so in novel situations, you have to look at the source (or .f90 files) to see what's what in your case.

[nobuhitomori]

Now I found my mistake. The atm pressure term

Code: Select all

P(i,j,N(ng)) = P(i,j,N(ng)) + 100._r8*Pair(i,j)/rho0

was correctly worked but my open boundary conditions
for storm surge condition got mass flux abnormally.
I checked above scheme qualitatively and it works
pretty fine.

I appreciate kate and other people's valuable comments.

[drews]

I'm using the 2-D ROMS configuration to model typhoon-driven storm surge. My model is forced with winds and atmospheric pressure. The code edits provided by rongzr produce a 1.5-meter rise in sea level for 850 mb, which is what I expect (after some sloshing). Thanks!

Here are the headers for my NetCDF forcing file, in case somebody else is trying to create their own:

netcdf winds { // and pressure
dimensions:
eta_u = 371 ;
xi_u = 250 ;
eta_v = 370 ;
xi_v = 251 ;
ocean_time = UNLIMITED ; // (3 currently)
eta_rho = 371 ;
xi_rho = 251 ;
variables:
double ocean_time(ocean_time) ;
ocean_time:long_name = "time since initialization" ;
ocean_time:units = "seconds since 0001-01-01 00:00:00" ;
ocean_time:calendar = "julian" ;
ocean_time:field = "time, scalar, series" ;
float sustr(ocean_time, eta_u, xi_u) ;
sustr:long_name = "surface u-momentum stress" ;
sustr:units = "newton meter-2" ;
sustr:time = "ocean_time" ;
sustr:field = "surface u-momentum stress, scalar, series" ;
float svstr(ocean_time, eta_v, xi_v) ;
svstr:long_name = "surface v-momentum stress" ;
svstr:units = "newton meter-2" ;
svstr:time = "ocean_time" ;
svstr:field = "surface v-momentum stress, scalar, series" ;
float Pair(ocean_time, eta_rho, xi_rho) ;
Pair:long_name = "surface air pressure" ;
Pair:units = "millibar" ;
Pair:time = "ocean_time" ;
Pair:field = "atmospheric pressure, scalar, series" ;
}

It looks like the new term provided by rongzr is:

1/2 * om_v * 100/rho0 * delta-depth * delta-Pair

Can somebody walk me through the physics of that term? What are on_u and om_v?

Carl

[rongzr]

You may want to check the 2D equation of ROMS. In the 2D case, adding atmos. pressure is simply add it to the pressure gradient term.