Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

STY: Use mbar and Pa for units of pressure #356

Open
wants to merge 5 commits into
base: dev/emc
Choose a base branch
from
Open

STY: Use mbar and Pa for units of pressure #356

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
50d209d
STY: Make capitalization of Pa consistent.
DWesl Jun 28, 2024
2965874
STY: Write millibar as mbar not mb.
DWesl Jun 28, 2024
d05140c
Replace pa with Pa in unit strings in nggps_diag
DWesl Jun 4, 2024
df6e7b8
STY: Use Pa instead of pa as pressure units
DWesl Aug 9, 2024
cb2aec9
Merge branch 'dev/emc' of NOAA-GFDL/GFDL_atmos_cubed_sphere into sty-…
DWesl Aug 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 10 additions & 10 deletions driver/fvGFS/fv_nggps_diag.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -255,14 +255,14 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)

if( Atm(n)%flagstruct%hydrostatic ) then
id_pfhy = register_diag_field ( trim(file_name), 'pfhy', axes(1:3), Time, &
'hydrostatic pressure', 'pa', missing_value=missing_value )
'hydrostatic pressure', 'Pa', missing_value=missing_value )
if (id_pfhy>0) then
kstt_pfhy = nlevs+1; kend_pfhy = nlevs+npzo
nlevs = nlevs + npzo
endif
else
id_pfnh = register_diag_field ( trim(file_name), 'pfnh', axes(1:3), Time, &
'non-hydrostatic pressure', 'pa', missing_value=missing_value )
'non-hydrostatic pressure', 'Pa', missing_value=missing_value )
if (id_pfnh>0) then
kstt_pfnh = nlevs+1; kend_pfnh = nlevs+npzo
nlevs = nlevs + npzo
Expand All @@ -282,7 +282,7 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
endif

id_omga = register_diag_field ( trim(file_name), 'omga', axes(1:3), Time, &
'Vertical pressure velocity', 'pa/sec', missing_value=missing_value )
'Vertical pressure velocity', 'Pa/sec', missing_value=missing_value )
if (id_omga>0) then
kstt_omga = nlevs+1; kend_omga = nlevs+npzo
nlevs = nlevs + npzo
Expand All @@ -296,7 +296,7 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
endif

id_delp = register_diag_field ( trim(file_name), 'delp', axes(1:3), Time, &
'pressure thickness', 'pa', missing_value=missing_value )
'pressure thickness', 'Pa', missing_value=missing_value )
if (id_delp>0) then
kstt_delp = nlevs+1; kend_delp = nlevs+npzo
nlevs = nlevs + npzo
Expand Down Expand Up @@ -325,7 +325,7 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
enddo
!
id_ps = register_diag_field ( trim(file_name), 'ps', axes(1:2), Time, &
'surface pressure', 'pa', missing_value=missing_value )
'surface pressure', 'Pa', missing_value=missing_value )
if( id_ps > 0) then
kstt_ps = nlevs+1; kend_ps = nlevs+1
nlevs = nlevs + 1
Expand Down Expand Up @@ -1374,7 +1374,7 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
endif
if( id_pfnh>0 ) then
call find_outputname(trim(file_name),'pfnh',output_name)
call add_field_to_bundle(trim(output_name),'non-hydrostatic pressure', 'pa', "time: point", &
call add_field_to_bundle(trim(output_name),'non-hydrostatic pressure', 'Pa', "time: point", &
axes(1:3), fcst_grid, kstt_pfnh,kend_pfnh, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
Expand All @@ -1387,15 +1387,15 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
else
if( id_pfhy>0 ) then
call find_outputname(trim(file_name),'pfhy',output_name)
call add_field_to_bundle(trim(output_name),'hydrostatic pressure', 'pa', "time: point", &
call add_field_to_bundle(trim(output_name),'hydrostatic pressure', 'Pa', "time: point", &
axes(1:3), fcst_grid, kstt_pfhy,kend_pfhy, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
endif
!
if( id_omga>0 ) then
call find_outputname(trim(file_name),'omga',output_name)
call add_field_to_bundle(trim(output_name),'Vertical pressure velocity', 'pa/sec', "time: point", &
call add_field_to_bundle(trim(output_name),'Vertical pressure velocity', 'Pa/sec', "time: point", &
axes(1:3), fcst_grid, kstt_omga,kend_omga, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
Expand All @@ -1410,7 +1410,7 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
!
if( id_delp > 0) then
call find_outputname(trim(file_name),'delp',output_name)
call add_field_to_bundle(trim(output_name),'pressure thickness', 'pa', "time: point", &
call add_field_to_bundle(trim(output_name),'pressure thickness', 'Pa', "time: point", &
axes(1:3), fcst_grid, kstt_delp,kend_delp, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
Expand Down Expand Up @@ -1458,7 +1458,7 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
!
if( id_ps > 0) then
call find_outputname(trim(file_name),'ps',output_name)
call add_field_to_bundle(trim(output_name),'surface pressure', 'pa', "time: point", &
call add_field_to_bundle(trim(output_name),'surface pressure', 'Pa', "time: point", &
axes(1:2), fcst_grid, kstt_ps,kend_ps, dyn_bundle, output_file, rcd=rc)
if(rc==0) num_field_dyn=num_field_dyn+1
endif
Expand Down
2 changes: 1 addition & 1 deletion tools/external_ic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3259,7 +3259,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
endif
endif ! data source /= FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE

! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
! For GFS spectral input, omega in Pa/sec is stored as w in the input data so actual w(m/s) is calculated
! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
! For FV3GFS NEMSIO input, w is already in m/s (but the code reads in as omga) and just needs to be remapped
!-------------------------------------------------------------
Expand Down
4 changes: 2 additions & 2 deletions tools/fv_diag_column.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -350,7 +350,7 @@ subroutine debug_column(pt, delp, delz, u, v, w, q, npz, ncnst, sphum, nwat, zvi
write(unit, *)

write(unit,500) 'k', 'T', 'delp', 'delz', 'u', 'v', 'w', 'sphum', 'cond', 'pres', 'NHprime'!, 'pdry', 'NHpdry'
write(unit,500) ' ', 'K', 'mb', 'm', 'm/s', 'm/s', 'm/s', 'g/kg', 'g/kg', 'mb', 'mb'!, ! 'mb', 'mb'
write(unit,500) ' ', 'K', 'mbar', 'm', 'm/s', 'm/s', 'm/s', 'g/kg', 'g/kg', 'mbar', 'mbar'!, ! 'mbar', 'mbar'
500 format(A4, A7, A8, A6, A8, A8, A8, A8, A9, A9, A9)
if (hydrostatic) then
call mpp_error(NOTE, 'Hydrostatic debug sounding not yet supported')
Expand Down Expand Up @@ -432,7 +432,7 @@ subroutine debug_column_dyn(pt, delp, delz, u, v, w, q, heat_source, cappa, akap
write(unit, *)

write(unit,500) 'k', 'T', 'delp', 'delz', 'u', 'v', 'w', 'sphum', 'cond', 'pres', 'NHprime', 'heat'
write(unit,500) ' ', 'K', 'mb', 'm', 'm/s', 'm/s', 'm/s', 'g/kg', 'g/kg', 'mb', 'mb', 'K'
write(unit,500) ' ', 'K', 'mbar', 'm', 'm/s', 'm/s', 'm/s', 'g/kg', 'g/kg', 'mbar', 'mbar', 'K'
500 format(A4, A7, A8, A6, A8, A8, A8, A8, A9, A9, A9, A8)
if (hydrostatic) then
call mpp_error(NOTE, 'Hydrostatic debug sounding not yet supported')
Expand Down
32 changes: 16 additions & 16 deletions tools/fv_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -433,9 +433,9 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
deallocate(grid_xt, grid_yt)
deallocate(grid_x, grid_y )

id_phalf = diag_axis_init('phalf', phalf, 'mb', 'z', &
id_phalf = diag_axis_init('phalf', phalf, 'mbar', 'z', &
'ref half pressure level', direction=-1, set_name="dynamics")
id_pfull = diag_axis_init('pfull', pfull, 'mb', 'z', &
id_pfull = diag_axis_init('pfull', pfull, 'mbar', 'z', &
'ref full pressure level', direction=-1, set_name="dynamics", edges=id_phalf)

!---- register static fields -------
Expand Down Expand Up @@ -547,16 +547,16 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
enddo
end if

id_plev = diag_axis_init('plev', levs(1:nplev)*1.0, 'mb', 'z', &
id_plev = diag_axis_init('plev', levs(1:nplev)*1.0, 'mbar', 'z', &
'actual pressure level', direction=-1, set_name="dynamics")

axe2(1) = id_xt
axe2(2) = id_yt
axe2(3) = id_plev

id_plev_ave_edges = diag_axis_init('plev_ave_edges', levs_ave(1:nplev_ave+1)*1.0, 'mb', 'z', &
id_plev_ave_edges = diag_axis_init('plev_ave_edges', levs_ave(1:nplev_ave+1)*1.0, 'mbar', 'z', &
'averaging layer pressure interface', direction=-1, set_name="dynamics")
id_plev_ave = diag_axis_init('plev_ave', (levs_ave(1:nplev_ave)+levs_ave(2:nplev_ave+1))*0.5, 'mb', 'z', &
id_plev_ave = diag_axis_init('plev_ave', (levs_ave(1:nplev_ave)+levs_ave(2:nplev_ave+1))*0.5, 'mbar', 'z', &
'averaging layer pressure level', direction=-1, set_name="dynamics", edges=id_plev_ave_edges)

axe_ave(1) = id_xt
Expand Down Expand Up @@ -994,32 +994,32 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
! Sea-level-pressure
!-------------------
id_slp = register_diag_field (trim(field), 'slp', axes(1:2), Time, &
'sea-level pressure', 'mb', missing_value=missing_value, &
'sea-level pressure', 'mbar', missing_value=missing_value, &
range=slprange )
!----------------------------------
! Bottom level pressure for masking
!----------------------------------
id_pmask = register_diag_field (trim(field), 'pmask', axes(1:2), Time, &
'masking pressure at lowest level', 'mb', &
'masking pressure at lowest level', 'mbar', &
missing_value=missing_value )
!------------------------------------------
! Fix for Bottom level pressure for masking
!------------------------------------------
id_pmaskv2 = register_diag_field(TRIM(field), 'pmaskv2', axes(1:2), Time,&
& 'masking pressure at lowest level', 'mb', missing_value=missing_value)
& 'masking pressure at lowest level', 'mbar', missing_value=missing_value)

!-------------------
! Hurricane scales:
!-------------------
! Net effects: ~ intensity * freq
id_c15 = register_diag_field (trim(field), 'cat15', axes(1:2), Time, &
'de-pression < 1000', 'mb', missing_value=missing_value)
'de-pression < 1000', 'mbar', missing_value=missing_value)
id_c25 = register_diag_field (trim(field), 'cat25', axes(1:2), Time, &
'de-pression < 980', 'mb', missing_value=missing_value)
'de-pression < 980', 'mbar', missing_value=missing_value)
id_c35 = register_diag_field (trim(field), 'cat35', axes(1:2), Time, &
'de-pression < 964', 'mb', missing_value=missing_value)
'de-pression < 964', 'mbar', missing_value=missing_value)
id_c45 = register_diag_field (trim(field), 'cat45', axes(1:2), Time, &
'de-pression < 944', 'mb', missing_value=missing_value)
'de-pression < 944', 'mbar', missing_value=missing_value)
! Frequency:
id_f15 = register_diag_field (trim(field), 'f15', axes(1:2), Time, &
'Cat15 frequency', 'none', missing_value=missing_value)
Expand Down Expand Up @@ -1063,16 +1063,16 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
'Relative Humidity', '%', missing_value=missing_value )
! 'Relative Humidity', '%', missing_value=missing_value, range=rhrange )
id_delp = register_diag_field ( trim(field), 'delp', axes(1:3), Time, &
'pressure thickness', 'pa', missing_value=missing_value )
'pressure thickness', 'Pa', missing_value=missing_value )
if ( .not. Atm(n)%flagstruct%hydrostatic ) &
id_delz = register_diag_field ( trim(field), 'delz', axes(1:3), Time, &
'height thickness', 'm', missing_value=missing_value )
if( Atm(n)%flagstruct%hydrostatic ) then
id_pfhy = register_diag_field ( trim(field), 'pfhy', axes(1:3), Time, &
'hydrostatic pressure', 'pa', missing_value=missing_value )
'hydrostatic pressure', 'Pa', missing_value=missing_value )
else
id_pfnh = register_diag_field ( trim(field), 'pfnh', axes(1:3), Time, &
'non-hydrostatic pressure', 'pa', missing_value=missing_value )
'non-hydrostatic pressure', 'Pa', missing_value=missing_value )
endif
id_zratio = register_diag_field ( trim(field), 'zratio', axes(1:3), Time, &
'nonhydro_ratio', 'n/a', missing_value=missing_value )
Expand Down Expand Up @@ -4607,7 +4607,7 @@ subroutine get_pressure_given_height(is, ie, js, je, ng, km, wz, kd, height, &
real, intent(in):: ts(is-ng:ie+ng,js-ng:je+ng)
real, intent(in):: peln(is:ie,km+1,js:je)
real, intent(in):: height(kd) !< must be monotonically decreasing with increasing k
real, intent(out):: a2(is:ie,js:je,kd) !< pressure (pa)
real, intent(out):: a2(is:ie,js:je,kd) !< pressure (Pa)
real, optional, intent(in):: fac

! local:
Expand Down
16 changes: 8 additions & 8 deletions tools/fv_eta.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -963,7 +963,7 @@ subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate)
do k=ks+1,km
tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) )
enddo
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hPa)=', tmp1/100.
write(*,800) (pm(k), k=km,1,-1)
endif

Expand Down Expand Up @@ -1137,7 +1137,7 @@ subroutine var_gfs(km, ak, bk, ptop, ks, pint, s_rate)
do k=ks+1,km
tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) )
enddo
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hPa)=', tmp1/100.
endif

end subroutine var_gfs
Expand Down Expand Up @@ -1313,7 +1313,7 @@ subroutine var_hi(km, ak, bk, ptop, ks, pint, s_rate)
do k=ks+1,km
tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) )
enddo
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hPa)=', tmp1/100.
endif


Expand Down Expand Up @@ -1470,7 +1470,7 @@ subroutine var_hi2(km, ak, bk, ptop, ks, pint, s_rate)
do k=ks+1,km
tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) )
enddo
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hPa)=', tmp1/100.
endif


Expand Down Expand Up @@ -1631,7 +1631,7 @@ subroutine var_dz(km, ak, bk, ptop, ks, pint, s_rate)
do k=ks+1,km
tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) )
enddo
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hPa)=', tmp1/100.
endif


Expand Down Expand Up @@ -1794,7 +1794,7 @@ subroutine var55_dz(km, ak, bk, ptop, ks, pint, s_rate)
do k=ks+1,km
tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) )
enddo
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hPa)=', tmp1/100.
endif


Expand Down Expand Up @@ -2366,7 +2366,7 @@ subroutine gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
enddo

ptop = pe1(1)
! if ( is_master() ) write(*,*) 'GW_1D: computed model top (pa)=', ptop
! if ( is_master() ) write(*,*) 'GW_1D: computed model top (Pa)=', ptop

! Set up "sigma" coordinate
ak(1) = pe1(1)
Expand Down Expand Up @@ -2519,7 +2519,7 @@ subroutine mount_waves(km, ak, bk, ptop, ks, pint)
do k=ks+1,km
tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) )
enddo
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hPa)=', tmp1/100.
endif

end subroutine mount_waves
Expand Down
14 changes: 7 additions & 7 deletions tools/fv_nudge.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,7 @@ module fv_nwp_nudge_mod
integer :: k_breed = 0
integer :: k_trop = 0
real :: p_trop = 950.E2
real :: dps_min = 50. !< maximum PS increment (pa; each step) due to inline breeding
real :: dps_min = 50. !< maximum PS increment (Pa; each step) due to inline breeding
real :: del2_cd = 0.16

real, allocatable:: s2c(:,:,:)
Expand Down Expand Up @@ -258,8 +258,8 @@ module fv_nwp_nudge_mod
real :: tau_vt_rad = 4.0

real :: pt_lim = 0.2
real :: slp_env = 101010. !< storm environment pressure (pa)
real :: pre0_env = 100000. !< critical storm environment pressure (pa) for size computation
real :: slp_env = 101010. !< storm environment pressure (Pa)
real :: pre0_env = 100000. !< critical storm environment pressure (Pa) for size computation
real, parameter:: tm_max = 315.
!------------------
real:: r_lo = 2.0
Expand Down Expand Up @@ -2598,7 +2598,7 @@ subroutine breed_slp_inline(nstep, dt, npz, ak, bk, phis, pe, pk, peln, pkz, del

call mp_reduce_sum(p_sum)
mass_sink = mass_sink / p_sum ! mean delta pressure to be added back to the environment to conserve mass
if(master .and. nudge_debug) write(*,*) 'TC#',n, 'Mass tele-ported (pa)=', mass_sink
if(master .and. nudge_debug) write(*,*) 'TC#',n, 'Mass tele-ported (Pa)=', mass_sink

!$OMP parallel do default(none) shared(is,ie,js,je,dist,r3,r2,ak,k_breed,delp,ps,mass_sink,npz) &
!$OMP private(pbreed, f1)
Expand Down Expand Up @@ -3269,8 +3269,8 @@ subroutine get_slp_obs(time, nobs, lon_obs, lat_obs, w10, mslp, slp_out, r_out,
integer, intent(in):: nobs !< number of observations in this particular storm
real(KIND=4), intent(in):: lon_obs(nobs)
real(KIND=4), intent(in):: lat_obs(nobs)
real(KIND=4), intent(in):: w10(nobs) !< observed 10-m widn speed
real(KIND=4), intent(in):: mslp(nobs) !< observed SLP in pa
real(KIND=4), intent(in):: w10(nobs) !< observed 10-m wind speed
real(KIND=4), intent(in):: mslp(nobs) !< observed SLP in Pa
real(KIND=4), intent(in):: slp_out(nobs) !< slp at r_out
real(KIND=4), intent(in):: r_out(nobs)
real(KIND=4), intent(in):: time_obs(nobs)
Expand All @@ -3279,7 +3279,7 @@ subroutine get_slp_obs(time, nobs, lon_obs, lat_obs, w10, mslp, slp_out, r_out,
! Output
real(kind=R_GRID), intent(out):: x_o , y_o !< position of the storm center
real, intent(out):: w10_o !< 10-m wind speed
real, intent(out):: slp_o !< Observed sea-level-pressure (pa)
real, intent(out):: slp_o !< Observed sea-level-pressure (Pa)
real, intent(out):: r_vor, p_vor
! Internal:
real:: t_thresh
Expand Down
4 changes: 2 additions & 2 deletions tools/test_cases.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2177,7 +2177,7 @@ subroutine init_case(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak,
enddo
ze1(1) = ztop

if ( is_master() ) write(*,*) 'Model top (pa)=', ptop
if ( is_master() ) write(*,*) 'Model top (Pa)=', ptop

do j=jsd,jed
do i=isd,ied
Expand Down Expand Up @@ -2246,7 +2246,7 @@ subroutine init_case(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak,
enddo
ze1(1) = ztop

if ( is_master() ) write(*,*) 'Model top (pa)=', ptop
if ( is_master() ) write(*,*) 'Model top (Pa)=', ptop

do j=jsd,jed
do i=isd,ied
Expand Down