apply the atmosphere pressure in Storm simulation
apply the atmosphere pressure in Storm simulation
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.
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.
I have been merging changes from Paul Budgell's code and noticed that he is doing this. He changed prsgrd32.h to include:
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.#ifdef SLP_GRAD
P(i,j,N(ng)) = P(i,j,N(ng)) + 100._r8*Pair(i,j)/rho0
#endif

 Posts: 22
 Joined: Fri Jul 08, 2005 5:42 pm
 Location: Kyoto University
Re: apply the atmosphere pressure in Storm simulation
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
as
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?
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)))
Code: Select all
P(i,j,N(ng)) = P(i,j,N(ng)) + 100._r8*Pair(i,j)/rho0
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?
Re: apply the atmosphere pressure in Storm simulation
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.
For the other, ROMS wants to read Pair in millibars. This looks like it would convert Pair to Pascals for this calculation.
Re: apply the atmosphere pressure in Storm simulation
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
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

 Posts: 22
 Joined: Fri Jul 08, 2005 5:42 pm
 Location: Kyoto University
Re: apply the atmosphere pressure in Storm simulation
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.
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.
Re: apply the atmosphere pressure in Storm simulation
According to my experience, you need to edited this two files
(1)mode_forces.FDefine 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 righthandside 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(i1,j)+ &
& h(i ,j))* &
& (gzeta(i1,j) &
& gzeta(i ,j))+ &
# if defined VAR_RHO_2D && defined SOLVE3D
& (h(i1,j) &
& h(i ,j))* &
& (gzetaSA(i1,j)+ &
& gzetaSA(i ,j)+ &
& cff2*(rhoA(i1,j) &
& rhoA(i ,j))* &
& (zwrk(i1,j) &
& zwrk(i ,j)))+ &
# endif
& (gzeta2(i1,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(i1,j)+h(i,j)+ &
& gzeta(i1,j)+gzeta(i,j))* &
& (Pair(i1,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,j1)+ &
& h(i,j ))* &
& (gzeta(i,j1) &
& gzeta(i,j ))+ &
# if defined VAR_RHO_2D && defined SOLVE3D
& (h(i,j1) &
& h(i,j ))* &
& (gzetaSA(i,j1)+ &
& gzetaSA(i,j )+ &
& cff2*(rhoA(i,j1) &
& rhoA(i,j ))* &
& (zwrk(i,j1) &
& zwrk(i,j )))+ &
# endif
& (gzeta2(i,j1) &
& 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,j1)+h(i,j)+ &
& gzeta(i,j1)+gzeta(i,j))* &
& (Pair(i,j1)Pair(i,j)) + &
# endif
# endif
& 0.0_r8
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note that this only works for the 2D case.
(1)mode_forces.FDefine 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 righthandside 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(i1,j)+ &
& h(i ,j))* &
& (gzeta(i1,j) &
& gzeta(i ,j))+ &
# if defined VAR_RHO_2D && defined SOLVE3D
& (h(i1,j) &
& h(i ,j))* &
& (gzetaSA(i1,j)+ &
& gzetaSA(i ,j)+ &
& cff2*(rhoA(i1,j) &
& rhoA(i ,j))* &
& (zwrk(i1,j) &
& zwrk(i ,j)))+ &
# endif
& (gzeta2(i1,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(i1,j)+h(i,j)+ &
& gzeta(i1,j)+gzeta(i,j))* &
& (Pair(i1,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,j1)+ &
& h(i,j ))* &
& (gzeta(i,j1) &
& gzeta(i,j ))+ &
# if defined VAR_RHO_2D && defined SOLVE3D
& (h(i,j1) &
& h(i,j ))* &
& (gzetaSA(i,j1)+ &
& gzetaSA(i,j )+ &
& cff2*(rhoA(i,j1) &
& rhoA(i,j ))* &
& (zwrk(i,j1) &
& zwrk(i,j )))+ &
# endif
& (gzeta2(i,j1) &
& 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,j1)+h(i,j)+ &
& gzeta(i,j1)+gzeta(i,j))* &
& (Pair(i,j1)Pair(i,j)) + &
# endif
# endif
& 0.0_r8
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note that this only works for the 2D case.
Re: apply the atmosphere pressure in Storm simulation
It's true that Paul's method only works in 3D. 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.

 Posts: 22
 Joined: Fri Jul 08, 2005 5:42 pm
 Location: Kyoto University
Re: apply the atmosphere pressure in Storm simulation
Now I found my mistake. The atm pressure term
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.
Code: Select all
P(i,j,N(ng)) = P(i,j,N(ng)) + 100._r8*Pair(i,j)/rho0
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
 Posts: 35
 Joined: Tue Jun 19, 2007 3:32 pm
 Location: National Center for Atmospheric Research
 Contact:
Re: apply the atmosphere pressure in Storm simulation
I'm using the 2D ROMS configuration to model typhoondriven storm surge. My model is forced with winds and atmospheric pressure. The code edits provided by rongzr produce a 1.5meter 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 00010101 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 umomentum stress" ;
sustr:units = "newton meter2" ;
sustr:time = "ocean_time" ;
sustr:field = "surface umomentum stress, scalar, series" ;
float svstr(ocean_time, eta_v, xi_v) ;
svstr:long_name = "surface vmomentum stress" ;
svstr:units = "newton meter2" ;
svstr:time = "ocean_time" ;
svstr:field = "surface vmomentum 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 * deltadepth * deltaPair
Can somebody walk me through the physics of that term? What are on_u and om_v?
Carl
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 00010101 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 umomentum stress" ;
sustr:units = "newton meter2" ;
sustr:time = "ocean_time" ;
sustr:field = "surface umomentum stress, scalar, series" ;
float svstr(ocean_time, eta_v, xi_v) ;
svstr:long_name = "surface vmomentum stress" ;
svstr:units = "newton meter2" ;
svstr:time = "ocean_time" ;
svstr:field = "surface vmomentum 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 * deltadepth * deltaPair
Can somebody walk me through the physics of that term? What are on_u and om_v?
Carl
Re: apply the atmosphere pressure in Storm simulation
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.