Doing some steps of the pressure calculation in double precision sinc…

…e the difference in output from single precision.

scm/src/scm_vgrid.F90

@@ -3,7 +3,7 @@
 
 module scm_vgrid
 
-use scm_kinds, only: sp, dp, qp
+use scm_kinds, only: sp, dp, qp, kind_scm_dp
 use scm_physical_constants, only : con_cp, con_rocp, con_fvirt, con_g, con_rd
 
 implicit none
@@ -48,6 +48,9 @@ subroutine get_FV3_vgrid(scm_input, scm_state)
 
       character(len=80)     :: line
       character(len=16)     :: file_format
+      integer :: nx, ny
+      real(kind_scm_dp), allocatable :: pres_l_row(:), pres_i(:,:)
+      real(kind_scm_dp), parameter :: zero = 0.0
 
 #include "fv_eta.h"
 
@@ -539,17 +542,24 @@ subroutine get_FV3_vgrid(scm_input, scm_state)
 
       p_ref = scm_input%input_pres_surf(1)
       pres_sfc_inv = 1.0/p_ref
+
+      nx = size(scm_state%pres_i, 1)
+      ny = size(scm_state%pres_i, 2)
+      allocate(pres_i(nx,ny), source=zero)
       do k=1, km+1
-        scm_state%pres_i(:,k) = scm_state%a_k(k) + scm_state%b_k(k)*p_ref
+        pres_i(:,k) = scm_state%a_k(k) + scm_state%b_k(k)*p_ref
         scm_state%si(:,k) = scm_state%a_k(k)*pres_sfc_inv + scm_state%b_k(k)
         scm_state%exner_i(:,k) = (scm_state%pres_i(:,k)/1.0E5)**con_rocp
       end do
+      scm_state%pres_i = pres_i
 
       !> - Calculate layer center pressures, sigma, and exner function.
+      allocate(pres_l_row(nx), source=zero)
       do k=1, km
-        scm_state%pres_l(:,k) = ((1.0/(con_rocp+1.0))*&
-          (scm_state%pres_i(:,k)**(con_rocp+1.0) - scm_state%pres_i(:,k+1)**(con_rocp+1.0))/ &
-          (scm_state%pres_i(:,k) - scm_state%pres_i(:,k+1)))**(1.0/con_rocp)
+        pres_l_row = ((1.0/(con_rocp+1.0))*&
+          (pres_i(:,k)**(con_rocp+1.0) - pres_i(:,k+1)**(con_rocp+1.0))/ &
+          (pres_i(:,k) - pres_i(:,k+1)))**(1.0/con_rocp)
+        scm_state%pres_l(:,k) = pres_l_row
         scm_state%sl(:,k) = 0.5*(scm_state%si(:,k) + scm_state%si(:,k+1))
 
         scm_state%exner_l(:,k) = (scm_state%pres_l(:,k)/1.0E5)**con_rocp