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

* STY: Make capitalization of Pa consistent.

It should keep my unit software from thinking I measure pressure in picoannum

* STY: Write millibar as mbar not mb.

My unit software thinks the latter is millibarn, a unit of area for atomic scattering.
This may help with CF compliance.

* Replace pa with Pa in unit strings in nggps_diag

My unit software reads "pa" as "picoannum", a unit of time on the order of microseconds, rather than as "pascals", the unit of pressure.

* STY: Use Pa instead of pa as pressure units

driver/UFS/fv_nggps_diag.F90

@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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

tools/external_ic.F90

@@ -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
 !-------------------------------------------------------------

tools/fv_diag_column.F90

@@ -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')
@@ -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')

tools/fv_diagnostics.F90

@@ -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 -------
@@ -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
@@ -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)
@@ -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 )
@@ -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:

tools/fv_eta.F90

@@ -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
 
@@ -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
@@ -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
 
 
@@ -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
 
 
@@ -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
 
 
@@ -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
 
 
@@ -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)
@@ -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

tools/fv_nudge.F90

@@ -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(:,:,:)
@@ -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
@@ -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)
@@ -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)
@@ -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

tools/test_cases.F90

@@ -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
@@ -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