diff --git a/components/cam/cime_config/buildnml b/components/cam/cime_config/buildnml
index d7e651efdb9..76bd6d38ef4 100755
--- a/components/cam/cime_config/buildnml
+++ b/components/cam/cime_config/buildnml
@@ -59,6 +59,15 @@ chdir "$CASEBUILD/camconf";
if ($BUILD_COMPLETE eq 'FALSE') {
+ # The following translation is hard-wired for backwards compatibility
+ # to support the differences between how the scripts specify the land grid
+ # and how it is specified internally
+
+ if ($ATM_GRID eq 'T31'){$ATM_GRID = "48x96";}
+ if ($ATM_GRID eq 'T42'){$ATM_GRID = "64x128";}
+ if ($ATM_GRID eq 'T85'){$ATM_GRID = "128x256";}
+ if ($ATM_GRID eq 'T341'){$ATM_GRID = "512x1024";}
+
# Some settings for single column mode.
my $scm = "";
if ($PTS_MODE eq 'TRUE') {$scm = "-scam -nosmp";}
diff --git a/components/cam/cime_config/config_compsets.xml b/components/cam/cime_config/config_compsets.xml
index fc0275deed2..6294bd395c1 100644
--- a/components/cam/cime_config/config_compsets.xml
+++ b/components/cam/cime_config/config_compsets.xml
@@ -526,6 +526,11 @@
FG20TRCN
20TR_CAM4_CLM40%CN_CICE%PRES_DOCN%DOM_RTM_CISM1_SWAV
+
+
+ F_ARM97_SCAM5
+ AR97_CAM5%SCAM_CLM40%SP_CICE%PRES_DOCN%DOM_RTM_SGLC_SWAV
+
@@ -551,7 +556,7 @@
1950-01-01
2004-01-01
1995-07-18
- 1997-06-18
+ 1997-06-19
@@ -595,14 +600,22 @@
ndays
+ ndays
2
+ 29
+
+
+
+ USGS
+
+
diff --git a/components/cam/src/control/vrtmap.F90 b/components/cam/src/control/vrtmap.F90
index 002cdf7a44b..5097287e9c5 100644
--- a/components/cam/src/control/vrtmap.F90
+++ b/components/cam/src/control/vrtmap.F90
@@ -1,3 +1,12 @@
+module vrtmap_mod
+
+implicit none
+
+private
+
+public :: vrtmap
+
+contains
subroutine vrtmap (pkdim ,pmap ,sigln ,dsigln ,kdpmap )
!-----------------------------------------------------------------------
@@ -12,14 +21,9 @@ subroutine vrtmap (pkdim ,pmap ,sigln ,dsigln ,kdpmap )
! Author: Jerry Olson
!
!-----------------------------------------------------------------------
- use shr_kind_mod, only: r8 => shr_kind_r8
+ use shr_kind_mod, only: r8 => shr_kind_r8
use cam_abortutils, only: endrun
- use cam_logfile, only: iulog
-#if (!defined UNICOSMP)
- use srchutil, only: ismin
-#endif
-!-----------------------------------------------------------------------
- implicit none
+ use cam_logfile, only: iulog
!-----------------------------------------------------------------------
!
! Arguments
@@ -42,15 +46,12 @@ subroutine vrtmap (pkdim ,pmap ,sigln ,dsigln ,kdpmap )
real(r8) del ! artificial grid interval
real(r8) dp ! artificial departure point
real(r8) eps ! epsilon factor
-#if (defined UNICOSMP)
- integer, external :: ismin
-#endif
!
!-----------------------------------------------------------------------
!
eps = 1.e-05_r8
del = ( sigln(pkdim) - sigln(1) )/real(pmap,r8)
- imin = ismin( pkdim-1,dsigln, 1 )
+ imin = minloc( dsigln(:pkdim-1), dim=1 )
if (del + eps >= dsigln(imin)) then
newmap = ( sigln(pkdim) - sigln(1) )/dsigln(imin) + 1
write(iulog,9000) pmap,newmap
@@ -73,3 +74,4 @@ subroutine vrtmap (pkdim ,pmap ,sigln ,dsigln ,kdpmap )
' Reset parameter "pmap" to at least ',i20)
end subroutine vrtmap
+end module vrtmap_mod
diff --git a/components/cam/src/dynamics/eul/dyn_comp.F90 b/components/cam/src/dynamics/eul/dyn_comp.F90
index a8200938cc1..b80e90df9aa 100644
--- a/components/cam/src/dynamics/eul/dyn_comp.F90
+++ b/components/cam/src/dynamics/eul/dyn_comp.F90
@@ -11,7 +11,7 @@ module dyn_comp
use constituents, only: sflxnam, tendnam, fixcnam, tottnam, hadvnam, vadvnam, cnst_get_ind
use pmgrid, only: plev, plevp, dyndecomp_set
use hycoef, only: hycoef_init
-use cam_history, only: dyn_decomp, addfld, add_default
+use cam_history, only: addfld, add_default, horiz_only
use phys_control, only: phys_getopts
use eul_control_mod, only: dyn_eul_readnl, eul_nsplit
use cam_logfile, only: iulog
@@ -36,11 +36,24 @@ module dyn_comp
integer :: placeholder
end type dyn_export_t
+! Frontogenesis indices
+integer, public :: frontgf_idx = -1
+integer, public :: frontga_idx = -1
+
+! Index into physics buffer for zonal mean zonal wind.
+integer, public :: uzm_idx = -1
+
!#######################################################################
CONTAINS
!#######################################################################
subroutine dyn_init(file, nlfilename)
+ use phys_control, only : use_gw_front
+ use dyn_grid, only: define_cam_grids, initgrid
+ use scamMod, only: single_column
+ use physics_buffer, only : pbuf_add_field, dtype_r8
+ use ppgrid, only : pcols, pver
+
! ARGUMENTS:
type(file_desc_t), intent(in) :: file ! PIO file handle for initial or restart file
character(len=*), intent(in) :: nlfilename
@@ -48,9 +61,10 @@ subroutine dyn_init(file, nlfilename)
! Local workspace
integer m ! Index
integer :: ixcldice, ixcldliq ! constituent indices for cloud liquid and ice water.
+ logical :: history_amwg ! output for AMWG diagnostics
logical :: history_budget ! output tendencies and state variables for CAM4
- ! temperature, water vapor, cloud ice and cloud
- ! liquid budgets.
+ ! temperature, water vapor, cloud ice and cloud
+ ! liquid budgets.
integer :: history_budget_histfile_num ! output history file number for budget fields
!----------------------------------------------------------------------------
@@ -64,57 +78,86 @@ subroutine dyn_init(file, nlfilename)
call spmdinit_dyn()
#endif
+ !----------------------------------------------------------------------
+
+ if (use_gw_front) then
+ call pbuf_add_field("FRONTGF", "global", dtype_r8, (/pcols,pver/), &
+ frontgf_idx)
+ call pbuf_add_field("FRONTGA", "global", dtype_r8, (/pcols,pver/), &
+ frontga_idx)
+ end if
+
! Initialize hybrid coordinate arrays
call hycoef_init(file)
- call addfld ('ETADOT ','1/s ',plevp,'A','Vertical (eta) velocity',dyn_decomp)
- call addfld ('U&IC ','m/s ',plev, 'I','Zonal wind' ,dyn_decomp )
- call addfld ('V&IC ','m/s ',plev, 'I','Meridional wind' ,dyn_decomp )
- call add_default ('U&IC ',0, 'I')
- call add_default ('V&IC ',0, 'I')
+ ! Run initgrid (the old initcom) which sets up coordinates and weights
+ call initgrid()
- call addfld ('PS&IC ','Pa ',1, 'I','Surface pressure' ,dyn_decomp )
- call addfld ('T&IC ','K ',plev, 'I','Temperature' ,dyn_decomp )
- call add_default ('PS&IC ',0, 'I')
- call add_default ('T&IC ',0, 'I')
+ ! Define the CAM grids (must be before addfld calls)
+ call define_cam_grids()
+
+ call addfld ('ETADOT',(/ 'ilev' /),'A', '1/s','Vertical (eta) velocity', gridname='gauss_grid')
+ call addfld ('U&IC', (/ 'lev' /), 'I', 'm/s','Zonal wind', gridname='gauss_grid' )
+ call addfld ('V&IC', (/ 'lev' /), 'I', 'm/s','Meridional wind', gridname='gauss_grid' )
+ call add_default ('U&IC',0, 'I')
+ call add_default ('V&IC',0, 'I')
+
+ call addfld ('PS&IC',horiz_only,'I', 'Pa','Surface pressure', gridname='gauss_grid' )
+ call addfld ('T&IC',(/ 'lev' /),'I', 'K','Temperature', gridname='gauss_grid' )
+ call add_default ('PS&IC',0, 'I')
+ call add_default ('T&IC',0, 'I')
do m = 1, pcnst
- call addfld (trim(cnst_name(m))//'&IC','kg/kg ',plev, 'I',cnst_longname(m) ,dyn_decomp )
+ call addfld (trim(cnst_name(m))//'&IC',(/ 'lev' /),'I', 'kg/kg',cnst_longname(m), gridname='gauss_grid' )
call add_default(trim(cnst_name(m))//'&IC',0, 'I')
- call addfld (hadvnam(m), 'kg/kg/s ',plev, 'A',trim(cnst_name(m))//' horizontal advection tendency ',dyn_decomp)
- call addfld (vadvnam(m), 'kg/kg/s ',plev, 'A',trim(cnst_name(m))//' vertical advection tendency ',dyn_decomp)
- call addfld (tendnam(m), 'kg/kg/s ',plev, 'A',trim(cnst_name(m))//' total tendency ',dyn_decomp)
- call addfld (tottnam(m), 'kg/kg/s ',plev, 'A',trim(cnst_name(m))//' horz + vert + fixer tendency ',dyn_decomp)
- call addfld (fixcnam(m), 'kg/kg/s ',plev, 'A',trim(cnst_name(m))//' tendency due to slt fixer',dyn_decomp)
+ call addfld (hadvnam(m), (/ 'lev' /), 'A', 'kg/kg/s',trim(cnst_name(m))//' horizontal advection tendency', &
+ gridname='gauss_grid')
+ call addfld (vadvnam(m), (/ 'lev' /), 'A', 'kg/kg/s',trim(cnst_name(m))//' vertical advection tendency', &
+ gridname='gauss_grid')
+ call addfld (tendnam(m), (/ 'lev' /), 'A', 'kg/kg/s',trim(cnst_name(m))//' total tendency', &
+ gridname='gauss_grid')
+ call addfld (tottnam(m), (/ 'lev' /), 'A', 'kg/kg/s',trim(cnst_name(m))//' horz + vert + fixer tendency', &
+ gridname='gauss_grid')
+ call addfld (fixcnam(m), (/ 'lev' /), 'A', 'kg/kg/s',trim(cnst_name(m))//' tendency due to slt fixer', &
+ gridname='gauss_grid')
end do
- call addfld ('DUH ','K/s ',plev, 'A','U horizontal diffusive heating',dyn_decomp)
- call addfld ('DVH ','K/s ',plev, 'A','V horizontal diffusive heating',dyn_decomp)
- call addfld ('DTH ','K/s ',plev, 'A','T horizontal diffusive heating',dyn_decomp)
+ call addfld ('DUH ',(/ 'lev' /),'A', 'K/s ','U horizontal diffusive heating', gridname='gauss_grid')
+ call addfld ('DVH ',(/ 'lev' /),'A', 'K/s ','V horizontal diffusive heating', gridname='gauss_grid')
+ call addfld ('DTH ',(/ 'lev' /),'A', 'K/s ','T horizontal diffusive heating', gridname='gauss_grid')
+
+ call addfld ('ENGYCORR',(/ 'lev' /),'A', 'W/m2 ','Energy correction for over-all conservation', gridname='gauss_grid')
+ call addfld ('TFIX ',horiz_only ,'A', 'K/s ','T fixer (T equivalent of Energy correction)', gridname='gauss_grid')
- call addfld ('ENGYCORR','W/m2 ',plev, 'A','Energy correction for over-all conservation',dyn_decomp)
- call addfld ('TFIX ','K/s ',1, 'A','T fixer (T equivalent of Energy correction)',dyn_decomp)
+ call addfld ('FU ',(/ 'lev' /),'A', 'm/s2 ','Zonal wind forcing term', gridname='gauss_grid')
+ call addfld ('FV ',(/ 'lev' /),'A', 'm/s2 ','Meridional wind forcing term', gridname='gauss_grid')
+ call addfld ('UTEND ',(/ 'lev' /),'A', 'm/s2 ','U tendency', gridname='gauss_grid')
+ call addfld ('VTEND ',(/ 'lev' /),'A', 'm/s2 ','V tendency', gridname='gauss_grid')
+ call addfld ('TTEND ',(/ 'lev' /),'A', 'K/s ','T tendency', gridname='gauss_grid')
+ call addfld ('LPSTEN ',horiz_only ,'A', 'Pa/s ','Surface pressure tendency', gridname='gauss_grid')
+ call addfld ('VAT ',(/ 'lev' /),'A', 'K/s ','Vertical advective tendency of T', gridname='gauss_grid')
+ call addfld ('KTOOP ',(/ 'lev' /),'A', 'K/s ','(Kappa*T)*(omega/P)', gridname='gauss_grid')
- call addfld ('FU ','m/s2 ',plev, 'A','Zonal wind forcing term',dyn_decomp)
- call addfld ('FV ','m/s2 ',plev, 'A','Meridional wind forcing term',dyn_decomp)
- call addfld ('UTEND ','m/s2 ',plev, 'A','U tendency',dyn_decomp)
- call addfld ('VTEND ','m/s2 ',plev, 'A','V tendency',dyn_decomp)
- call addfld ('TTEND ','K/s ',plev, 'A','T tendency',dyn_decomp)
- call addfld ('LPSTEN ','Pa/s ',1, 'A','Surface pressure tendency',dyn_decomp)
- call addfld ('VAT ','K/s ',plev, 'A','Vertical advective tendency of T',dyn_decomp)
- call addfld ('KTOOP ','K/s ',plev, 'A','(Kappa*T)*(omega/P)',dyn_decomp)
+ call phys_getopts(history_amwg_out=history_amwg, &
+ history_budget_out = history_budget, &
+ history_budget_histfile_num_out = history_budget_histfile_num)
+
+ if (history_amwg) then
+ call add_default ('DTH ', 1, ' ')
+ end if
- call add_default ('DTH ', 1, ' ')
- call phys_getopts(history_budget_out = history_budget, history_budget_histfile_num_out = history_budget_histfile_num)
if ( history_budget ) then
call cnst_get_ind('CLDLIQ', ixcldliq)
call cnst_get_ind('CLDICE', ixcldice)
- call add_default(hadvnam( 1), history_budget_histfile_num, ' ')
- call add_default(hadvnam(ixcldliq), history_budget_histfile_num, ' ')
- call add_default(hadvnam(ixcldice), history_budget_histfile_num, ' ')
- call add_default(vadvnam( 1), history_budget_histfile_num, ' ')
- call add_default(vadvnam(ixcldliq), history_budget_histfile_num, ' ')
- call add_default(vadvnam(ixcldice), history_budget_histfile_num, ' ')
+ ! The following variables are not defined for single column
+ if (.not. single_column) then
+ call add_default(hadvnam( 1), history_budget_histfile_num, ' ')
+ call add_default(hadvnam(ixcldliq), history_budget_histfile_num, ' ')
+ call add_default(hadvnam(ixcldice), history_budget_histfile_num, ' ')
+ call add_default(vadvnam( 1), history_budget_histfile_num, ' ')
+ call add_default(vadvnam(ixcldliq), history_budget_histfile_num, ' ')
+ call add_default(vadvnam(ixcldice), history_budget_histfile_num, ' ')
+ end if
call add_default(fixcnam( 1), history_budget_histfile_num, ' ')
call add_default(fixcnam(ixcldliq), history_budget_histfile_num, ' ')
call add_default(fixcnam(ixcldice), history_budget_histfile_num, ' ')
@@ -124,11 +167,11 @@ subroutine dyn_init(file, nlfilename)
call add_default(tendnam( 1), history_budget_histfile_num, ' ')
call add_default(tendnam(ixcldliq), history_budget_histfile_num, ' ')
call add_default(tendnam(ixcldice), history_budget_histfile_num, ' ')
- call add_default('TTEND ' , history_budget_histfile_num, ' ')
- call add_default('TFIX ' , history_budget_histfile_num, ' ')
- call add_default('KTOOP ' , history_budget_histfile_num, ' ')
- call add_default('VAT ' , history_budget_histfile_num, ' ')
- call add_default('DTH ' , history_budget_histfile_num, ' ')
+ call add_default('TTEND', history_budget_histfile_num, ' ')
+ call add_default('TFIX', history_budget_histfile_num, ' ')
+ call add_default('KTOOP', history_budget_histfile_num, ' ')
+ call add_default('VAT', history_budget_histfile_num, ' ')
+ call add_default('DTH', history_budget_histfile_num, ' ')
end if
end subroutine dyn_init
diff --git a/components/cam/src/dynamics/eul/dyn_grid.F90 b/components/cam/src/dynamics/eul/dyn_grid.F90
index c23ecdc2683..bde40252f14 100644
--- a/components/cam/src/dynamics/eul/dyn_grid.F90
+++ b/components/cam/src/dynamics/eul/dyn_grid.F90
@@ -40,6 +40,9 @@ module dyn_grid
implicit none
save
+ ! The Eulerian dynamics grids
+ integer, parameter, public :: dyn_decomp = 101
+
integer, private :: ngcols_d = 0 ! number of dynamics columns
integer, parameter :: ptimelevels = 3 ! number of time levels in the dycore
@@ -47,6 +50,341 @@ module dyn_grid
contains
+ subroutine initgrid
+
+ !-----------------------------------------------------------------------
+ !
+ ! Purpose:
+ ! Initialize Model commons, including COMCON, COMHYB, COMMAP, COMSPE,
+ ! and COMTRCNM
+ !
+ ! Method:
+ !
+ ! Author:
+ ! Original version: CCM1
+ ! Standardized: L. Bath, Jun 1992
+ ! L. Buja, Feb 1996
+ !
+ !-----------------------------------------------------------------------
+ !
+ ! $Id$
+ ! $Author$
+ !
+ !-----------------------------------------------------------------------
+ use shr_kind_mod, only: r8 => shr_kind_r8
+ use pmgrid, only: plat, plev, plon, plevp
+ use spmd_utils, only: masterproc, iam
+ use pspect
+ use comspe
+ use rgrid, only: nlon, beglatpair, wnummax, nmmax, fullgrid
+ use scanslt, only: nlonex, platd, j1
+ use gauaw_mod, only: gauaw
+ use commap, only: sq, rsq, slat, w, cs, href, ecref, clat, clon, &
+ latdeg, londeg, xm
+ use physconst, only: rair, rearth, ra
+ use time_manager, only: get_step_size
+ use scamMod, only: scmlat,scmlon,single_column
+ use hycoef, only: hypd, hypm
+ use eul_control_mod, only: ifax, trig,eul_nsplit
+ use cam_logfile, only: iulog
+ !-----------------------------------------------------------------------
+ implicit none
+ !-----------------------------------------------------------------------
+
+ ! Local workspace
+ !
+ real(r8) zsi(plat) ! sine of latitudes
+ real(r8) zw(plat) ! Gaussian weights
+ real(r8) zra2 ! ra squared
+ real(r8) zalp(2*pspt) ! Legendre function array
+ real(r8) zdalp(2*pspt) ! Derivative array
+ real(r8) zslat ! sin of lat and cosine of colatitude
+
+ integer i ! longitude index
+ integer j ! Latitude index
+ integer k ! Level index
+ integer kk ! Level index
+ integer kkk ! Level index
+ integer m,lm,mr,lmr ! Indices for legendre array
+ integer n ! Index for legendre array
+ integer nkk ! Print control variables
+ integer ik1 ! Print index temporary variable
+ integer ik2 ! Print index temporary variable
+ integer itmp ! Dimension of polynomial arrays temporary.
+ integer iter ! Iteration index
+ real(r8) zdt ! Time step for settau
+
+ logical lprint ! Debug print flag
+ integer irow ! Latitude pair index
+ integer lat ! Latitude index
+
+ real(r8) xlat ! Latitude (radians)
+ real(r8) pi ! Mathematical pi (3.14...)
+ real(r8) dtime ! timestep size [seconds]
+ !
+ !-----------------------------------------------------------------------
+ !
+ lprint = masterproc .and. .FALSE.
+
+ dtime = get_step_size()
+ zdt = dtime/eul_nsplit
+ if (masterproc) write(iulog,*) 'zdt, dtime, eul_nsplit=',zdt,dtime,eul_nsplit
+
+ ! call hdinti (rearth ,dtime )
+ call hdinti (rearth ,zdt )
+
+ if ( .not. single_column ) then
+ !
+ ! NMAX dependent arrays
+ !
+ if (pmmax.gt.plon/2) then
+ call endrun ('INITGRID:mmax=ptrm+1 .gt. plon/2')
+ end if
+ endif
+
+ zra2 = ra*ra
+ do j=2,pnmax
+ sq(j) = j*(j-1)*zra2
+ rsq(j) = 1._r8/sq(j)
+ end do
+ sq(1) = 0._r8
+ rsq(1) = 0._r8
+ !
+ ! MMAX dependent arrays
+ !
+ do j=1,pmmax
+ xm(j) = j-1
+ end do
+ if ( .not. single_column ) then
+ !
+ ! Gaussian latitude dependent arrays
+ !
+ call gauaw(zsi ,zw ,plat )
+ do irow=1,plat/2
+ slat(irow) = zsi(irow)
+ w(irow) = zw(irow)
+ w(plat - irow + 1) = zw(irow)
+ cs(irow) = 1._r8 - zsi(irow)*zsi(irow)
+ xlat = asin(slat(irow))
+ clat(irow) = -xlat
+ clat(plat - irow + 1) = xlat
+ end do
+
+ do lat=1,plat
+ latdeg(lat) = clat(lat)*45._r8/atan(1._r8)
+ end do
+ endif
+ !
+ ! Integration matrices of hydrostatic equation(href) and conversion
+ ! term(a). href computed as in ccm0 but isothermal bottom ecref
+ ! calculated to conserve energy
+ !
+ do k=1,plev
+ do kk=1,plev
+ href(kk,k) = 0._r8
+ ecref(kk,k) = 0._r8
+ end do
+ end do
+ !
+ ! Mean atmosphere energy conversion term is consistent with continiuty
+ ! Eq. In ecref, 1st index = column; 2nd index = row of matrix.
+ ! Mean atmosphere energy conversion term is energy conserving
+ !
+ do k=1,plev
+ ecref(k,k) = 0.5_r8/hypm(k) * hypd(k)
+ do kk=1,k-1
+ ecref(kk,k) = 1._r8/hypm(k) * hypd(kk)
+ end do
+ end do
+ !
+ ! Reference hydrostatic integration matrix consistent with conversion
+ ! term for energy conservation. In href, 1st index = column;
+ ! 2nd index = row of matrix.
+ !
+ do k = 1,plev
+ do kk = k,plev
+ href(kk,k) = ecref(k,kk)*hypd(kk)/hypd(k)
+ end do
+ end do
+ !
+ ! Print statements
+ !
+ if (lprint) then
+ nkk = plev/13
+ if (mod(plev,13).ne.0) nkk = nkk + 1
+ write(iulog,*)' '
+ write(iulog,*)'INITGRID: Hydrostatic matrix href'
+ do kk=1,nkk
+ ik1 = 1 + (kk-1)*13
+ ik2 = min0( ik1+12, plev )
+ write(iulog,9920) (k,k=ik1,ik2)
+ do kkk=1,plev
+ write(iulog,9910) kkk,(href(kkk,k),k=ik1,ik2)
+ end do
+ end do
+ write(iulog,*)' '
+ write(iulog,*)'INITGRID: Thermodynamic matrix ecref'
+ do kk=1,nkk
+ ik1 = 1 + (kk-1)*13
+ ik2 = min0( ik1+12, plev )
+ write(iulog,9920) (k,k=ik1,ik2)
+ do kkk=1,plev
+ write(iulog,9910) kkk,(ecref(kkk,k),k=ik1,ik2)
+ end do
+ end do
+ end if
+ !
+ ! Multiply href by r
+ !
+ do k=1,plev
+ do kk=1,plev
+ href(kk,k) = href(kk,k)*rair
+ end do
+ end do
+ !
+ ! Compute truncation parameters
+ !
+ if (masterproc) then
+ write(iulog,9950) ptrm,ptrn,ptrk
+ end if
+ !
+ ! Determine whether full or reduced grid
+ !
+ fullgrid = .true.
+ do j=1,plat
+ if (masterproc) then
+ write(iulog,*)'nlon(',j,')=',nlon(j),' wnummax(',j,')=',wnummax(j)
+ end if
+ if (nlon(j).lt.plon) fullgrid = .false.
+ end do
+
+ if ( single_column ) then
+ do j=1,plat
+ slat(j) = 1.0_r8 * sin(4.0_r8*atan(1.0_r8)*scmlat/180._r8)
+ w(j) = 2.0_r8/plat
+ cs(j) = 10._r8 - slat(j)*slat(j)
+ ! itmp = 2*pspt - 1
+ ! call phcs (zalp ,zdalp ,itmp ,zslat )
+ ! call reordp(j ,itmp ,zalp ,zdalp )
+ end do
+
+ !
+ ! Latitude array (list of latitudes in radians)
+ !
+ xlat = asin(slat(1))
+ clat(1) = xlat
+
+ clat(1)=scmlat*atan(1._r8)/45._r8
+ latdeg(1) = clat(1)*45._r8/atan(1._r8)
+ clon(1,1) = 4.0_r8*atan(1._r8)*mod((scmlon+360._r8),360._r8)/180._r8
+ londeg(1,1) = mod((scmlon+360._r8),360._r8)
+ !
+ ! SCAM not yet able to handle reduced grid.
+ !
+ if (.not. fullgrid) then
+ call endrun ('INITGRID: SCAM not yet configured to handle reduced grid')
+ end if
+ else
+ !
+ ! Compute constants related to Legendre transforms
+ ! Compute and reorder ALP and DALP
+ !
+ allocate( alp (pspt,plat/2) )
+ allocate( dalp (pspt,plat/2) )
+ do j=1,plat/2
+ zslat = slat(j)
+ itmp = 2*pspt - 1
+ call phcs (zalp ,zdalp ,itmp ,zslat )
+ call reordp(j ,itmp ,zalp ,zdalp )
+ end do
+ !
+ ! Copy and save local ALP and DALP
+ !
+ allocate( lalp (lpspt,plat/2) )
+ allocate( ldalp (lpspt,plat/2) )
+ do j=1,plat/2
+ do lm=1,numm(iam)
+ m = locm(lm,iam)
+ mr = nstart(m)
+ lmr = lnstart(lm)
+ do n=1,nlen(m)
+ lalp(lmr+n,j) = alp(mr+n,j)
+ ldalp(lmr+n,j) = dalp(mr+n,j)
+ end do
+ end do
+ end do
+ !
+ ! Mirror latitudes south of south pole
+ !
+ lat = 1
+ do j=j1-2,1,-1
+ nlonex(j) = nlon(lat)
+ lat = lat + 1
+ end do
+ nlonex(j1-1) = nlon(1) ! south pole
+ !
+ ! Real latitudes
+ !
+ j = j1
+ do lat=1,plat
+ nlonex(j) = nlon(lat)
+ j = j + 1
+ end do
+ nlonex(j1+plat) = nlon(plat) ! north pole
+ !
+ ! Mirror latitudes north of north pole
+ !
+ lat = plat
+ do j=j1+plat+1,platd
+ nlonex(j) = nlon(lat)
+ lat = lat - 1
+ end do
+ !
+ ! Longitude array
+ !
+ pi = 4.0_r8*atan(1.0_r8)
+ do lat=1,plat
+ do i=1,nlon(lat)
+ londeg(i,lat) = (i-1)*360._r8/nlon(lat)
+ clon(i,lat) = (i-1)*2.0_r8*pi/nlon(lat)
+ end do
+ end do
+
+ do j=1,plat/2
+ nmmax(j) = wnummax(j) + 1
+ end do
+ do m=1,pmmax
+ do irow=1,plat/2
+ if (nmmax(irow) .ge. m) then
+ beglatpair(m) = irow
+ goto 10
+ end if
+ end do
+ call endrun ('INITGRID: Should not ever get here')
+10 continue
+ end do
+ !
+ ! Set up trigonometric tables for fft
+ !
+ do j=1,plat
+ call set99(trig(1,j),ifax(1,j),nlon(j))
+ end do
+ !
+ ! Set flag indicating dynamics grid is now defined.
+ ! NOTE: this ASSUMES initgrid is called after spmdinit. The setting of nlon done here completes
+ ! the definition of the dynamics grid.
+ !
+ endif
+
+ return
+
+9910 format( 1x,i3,13f9.5)
+9920 format(/, 13i9)
+9950 format(/,' Truncation Parameters',/,' NTRM = ',i4,/, &
+ ' NTRN = ',i4,/,' NTRK = ',i4,/)
+
+ end subroutine initgrid
+
subroutine get_block_ldof_d(nlev, ldof)
use pio, only : pio_offset_kind
integer, intent(in) :: nlev
@@ -369,7 +707,7 @@ end subroutine get_horiz_grid_dim_d
!========================================================================
!
subroutine get_horiz_grid_d(size,clat_d_out,clon_d_out,area_d_out, &
- wght_d_out)
+ wght_d_out,lat_d_out,lon_d_out)
!-----------------------------------------------------------------------
!
@@ -397,6 +735,8 @@ subroutine get_horiz_grid_d(size,clat_d_out,clon_d_out,area_d_out, &
! area
real(r8), intent(out), optional :: wght_d_out(size) ! column integration
! weight
+ real(r8), intent(out), optional :: lat_d_out(:) ! column degree latitudes
+ real(r8), intent(out), optional :: lon_d_out(:) ! column degree longitudes
!---------------------------Local workspace-----------------------------
!
integer i,j ! loop indices
@@ -499,6 +839,40 @@ subroutine get_horiz_grid_d(size,clat_d_out,clon_d_out,area_d_out, &
enddo
endif
+ if(present(lon_d_out)) then
+ if(size == ngcols_d) then
+ n = 0
+ do j = 1,plat
+ do i = 1,nlon(j)
+ n = n + 1
+ lon_d_out(n) = londeg(i,j)
+ end do
+ end do
+ else if(size == plon) then
+ lon_d_out(:) = londeg(:,1)
+ else
+ write(iulog,*)'GET_HORIZ_GRID_D: arrays not large enough (', &
+ size,' < ',ngcols_d,' ) '
+ call endrun
+ end if
+ end if
+ if(present(lat_d_out)) then
+ if(size == ngcols_d) then
+ n = 0
+ do j = 1,plat
+ do i = 1,nlon(j)
+ n = n + 1
+ lat_d_out(n) = latdeg(j)
+ end do
+ end do
+ else if(size == plat) then
+ lat_d_out(:) = latdeg(:)
+ else
+ write(iulog,*)'GET_HORIZ_GRID_D: arrays not large enough (', &
+ size,' < ',ngcols_d,' ) '
+ call endrun
+ end if
+ end if
!
return
end subroutine get_horiz_grid_d
@@ -755,4 +1129,104 @@ subroutine dyn_grid_get_elem_coords( latndx, rlon, rlat, cdex )
!#######################################################################
+subroutine define_cam_grids()
+ use pspect, only: ptrm, ptrn, ptrk
+ use pmgrid, only: beglat, endlat, plon, plat
+ use commap, only: londeg, latdeg, w
+ use cam_grid_support, only: horiz_coord_t, horiz_coord_create, iMap
+ use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register
+
+ ! Local variables
+ integer :: i, j, ind
+ integer(iMap), pointer :: grid_map(:,:)
+ integer(iMap) :: latmap(endlat - beglat + 1)
+ type(horiz_coord_t), pointer :: lat_coord
+ type(horiz_coord_t), pointer :: lon_coord
+ real(r8), pointer :: rattval(:)
+
+ nullify(grid_map)
+ nullify(lat_coord)
+ nullify(lon_coord)
+ nullify(rattval)
+
+ ! Dynamics Grid
+ ! Make grid and lat maps (need to do this because lat indices are distributed)
+ ! Note that for this dycore, some pes may be inactive
+ if(endlat >= beglat) then
+ allocate(grid_map(4, (plon * (endlat - beglat + 1))))
+ ind = 0
+ do i = beglat, endlat
+ do j = 1, plon
+ ind = ind + 1
+ grid_map(1, ind) = j
+ grid_map(2, ind) = i
+ grid_map(3, ind) = j
+ grid_map(4, ind) = i
+ end do
+ end do
+ ! Do we need a lat map?
+ if ((beglat /= 1) .or. (endlat /= plat)) then
+ do i = beglat, endlat
+ latmap(i - beglat + 1) = i
+ end do
+ end if
+ else
+ allocate(grid_map(4, 0))
+ end if
+
+ ! Create the lat coordinate
+ if ((beglat /= 1) .or. (endlat /= plat)) then
+ lat_coord => horiz_coord_create('lat', '', plat, 'latitude', &
+ 'degrees_north', beglat, endlat, latdeg(beglat:endlat), map=latmap)
+ else
+ lat_coord => horiz_coord_create('lat', '', plat, 'latitude', &
+ 'degrees_north', beglat, endlat, latdeg(beglat:endlat))
+ end if
+
+ ! Create the lon coordinate
+ lon_coord => horiz_coord_create('lon', '', plon, 'longitude', &
+ 'degrees_east', 1, plon, londeg(1:plon, 1))
+
+ call cam_grid_register('gauss_grid', dyn_decomp, lat_coord, lon_coord, &
+ grid_map, unstruct=.false.)
+
+ allocate(rattval(size(w)))
+ rattval = w
+ call cam_grid_attribute_register('gauss_grid', 'gw', 'gauss weights', 'lat', rattval)
+ nullify(rattval) ! belongs to attribute
+
+ ! Scalar variable 'attributes'
+ call cam_grid_attribute_register('gauss_grid', 'ntrm', &
+ 'spectral truncation parameter M', ptrm)
+ call cam_grid_attribute_register('gauss_grid', 'ntrn', &
+ 'spectral truncation parameter N', ptrn)
+ call cam_grid_attribute_register('gauss_grid', 'ntrk', &
+ 'spectral truncation parameter K', ptrk)
+ ! These belong to the grid now
+ nullify(grid_map)
+ nullify(lat_coord)
+ nullify(lon_coord)
+
+end subroutine define_cam_grids
+
+!#######################################################################
+
+subroutine physgrid_copy_attributes_d(gridname, grid_attribute_names)
+ use cam_grid_support, only: max_hcoordname_len
+
+ ! Dummy arguments
+ character(len=max_hcoordname_len), intent(out) :: gridname
+ character(len=max_hcoordname_len), pointer, intent(out) :: grid_attribute_names(:)
+
+ gridname = 'gauss_grid'
+ allocate(grid_attribute_names(4))
+ grid_attribute_names(1) = 'gw'
+ grid_attribute_names(2) = 'ntrm'
+ grid_attribute_names(3) = 'ntrn'
+ grid_attribute_names(4) = 'ntrk'
+
+end subroutine physgrid_copy_attributes_d
+
+!#######################################################################
+
end module dyn_grid
diff --git a/components/cam/src/dynamics/eul/inidat.F90 b/components/cam/src/dynamics/eul/inidat.F90
index bb05f0ce3f3..961c975b085 100644
--- a/components/cam/src/dynamics/eul/inidat.F90
+++ b/components/cam/src/dynamics/eul/inidat.F90
@@ -8,26 +8,25 @@ module inidat
!
! Author: J. Olson May 2004
!
-! Modified: A. Gettelman and C. Craig Nov 2010 - put micro/macro physics into separate routines
-!
!-----------------------------------------------------------------------
- use ppgrid, only: begchunk, endchunk, pcols
- use pmgrid, only: beglat, endlat, plon, plat, plev, plnlv
- use rgrid, only: nlon
- use prognostics, only : div, vort, t3, u3, v3, q3, n3, phis, dpsm, dpsl, ps
- use ncdio_atm, only: infld
- use shr_kind_mod, only: r8 => shr_kind_r8
+ use ppgrid, only: begchunk, endchunk, pcols
+ use pmgrid, only: beglat, endlat, plon, plat, plev, plnlv
+ use rgrid, only: nlon
+ use prognostics, only : div, vort, t3, u3, v3, q3, n3, phis, dpsm, dpsl, ps
+ use ncdio_atm, only: infld
+ use shr_kind_mod, only: r8 => shr_kind_r8
use cam_abortutils , only: endrun
- use phys_grid, only: get_ncols_p
+ use phys_grid, only: get_ncols_p
- use spmd_utils, only: masterproc, mpicom, mpir8
- use cam_control_mod, only: ideal_phys, aqua_planet, moist_physics, adiabatic
- use scamMod, only: single_column, use_camiop, have_u, have_v, &
- have_cldliq, have_cldice,loniop,latiop,scmlat,scmlon
- use cam_logfile, only: iulog
- use pio, only: file_desc_t, pio_noerr, pio_inq_varid, pio_get_att, &
+ use spmd_utils, only: masterproc, mpicom, mpir8
+ use cam_control_mod, only: ideal_phys, aqua_planet, moist_physics, adiabatic
+ use cam_initfiles, only: initial_file_get_id, topo_file_get_id
+ use scamMod, only: single_column, use_camiop, have_u, have_v, &
+ have_cldliq, have_cldice,loniop,latiop,scmlat,scmlon
+ use cam_logfile, only: iulog
+ use pio, only: file_desc_t, pio_noerr, pio_inq_varid, pio_get_att, &
pio_inq_attlen, pio_inq_dimid, pio_inq_dimlen, pio_get_var,var_desc_t, &
- pio_seterrorhandling, pio_bcast_error, pio_internal_error
+ pio_seterrorhandling, pio_bcast_error, pio_internal_error, pio_offset_kind
implicit none
PRIVATE
@@ -60,7 +59,7 @@ module inidat
!*********************************************************************
- subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
+ subroutine read_inidat(fh_ini, fh_topo, dyn_in)
!
!-----------------------------------------------------------------------
!
@@ -69,21 +68,21 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
!
!-----------------------------------------------------------------------
!
- use ppgrid, only: pverp
- use phys_grid, only: scatter_field_to_chunk, gather_chunk_to_field,clat_p,clon_p
+ use phys_grid, only: clat_p, clon_p
use constituents, only: pcnst, cnst_name, cnst_read_iv, cnst_get_ind
use commap, only: clat,clon
use iop, only: setiopupdate,readiopdata
use dyn_comp , only: dyn_import_t
- use physconst, only: pi
-
+ use physconst, only: pi
+ use cam_pio_utils, only: cam_pio_get_var
+
!
! Arguments
!
- type(file_desc_t), intent(inout) :: ncid_ini, ncid_topo
+ implicit none
+ type(file_desc_t),intent(inout) :: fh_ini, fh_topo
type(dyn_import_t) :: dyn_in ! not used in eul dycore
- integer, optional, intent(in) :: datatype ! in: model or analysis dataset
integer type
!
!---------------------------Local workspace-----------------------------
@@ -91,6 +90,8 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
integer i,c,m,n,lat ! indices
integer ncol
+! type(file_desc_t), pointer :: fh_ini, fh_topo
+
real(r8), pointer, dimension(:,:,:) :: convptr_2d
real(r8), pointer, dimension(:,:,:,:) :: convptr_3d
real(r8), pointer, dimension(:,:,:,:) :: cldptr
@@ -104,6 +105,8 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
integer londimid,dimlon,latdimid,dimlat,latvarid,lonvarid
integer strt(3),cnt(3)
+ character(len=3), parameter :: arraydims3(3) = (/ 'lon', 'lev', 'lat' /)
+ character(len=3), parameter :: arraydims2(2) = (/ 'lon', 'lat' /)
type(var_desc_t) :: varid
real(r8), allocatable :: tmp2d(:,:)
!
@@ -122,8 +125,10 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
! processors on the physics grid.
!
!-----------------------------------------------------------------------
-!
-!
+
+! fh_ini => initial_file_get_id()
+! fh_topo => topo_file_get_id()
+
!-------------------------------------
! Allocate memory for temporary arrays
!-------------------------------------
@@ -150,21 +155,18 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
fieldname = 'U'
- call infld(fieldname, ncid_ini, 'lon', 'lev', 'lat', 1, plon, 1, plev, 1, plat, &
- arr3d_a, readvar, grid_map='global')
+ call cam_pio_get_var(fieldname, fh_ini, arraydims3, arr3d_a, found=readvar)
if(.not. readvar) call endrun('dynamics/eul/inidat.F90')
fieldname = 'V'
- call infld(fieldname, ncid_ini, 'lon', 'lev', 'lat', 1, plon, 1, plev, 1, plat, &
- arr3d_b, readvar, grid_map='global')
+ call cam_pio_get_var(fieldname, fh_ini, arraydims3, arr3d_b, found=readvar)
if(.not. readvar) call endrun('dynamics/eul/inidat.F90')
call process_inidat('UV')
fieldname = 'T'
- call infld(fieldname, ncid_ini, 'lon', 'lev', 'lat', 1, plon, 1, plev, 1, plat, &
- t3_tmp, readvar, grid_map='global')
+ call cam_pio_get_var(fieldname, fh_ini, arraydims3, t3_tmp, found=readvar)
if(.not. readvar) call endrun('dynamics/eul/inidat.F90')
call process_inidat('T')
@@ -175,10 +177,10 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
readvar = .false.
fieldname = cnst_name(m)
- if(cnst_read_iv(m)) &
- call infld(fieldname, ncid_ini, 'lon', 'lev', 'lat', 1, plon, 1, plev, 1, plat, &
- arr3d_a, readvar, grid_map='global')
- call process_inidat('CONSTS', m_cnst=m, ncid=ncid_ini)
+ if(cnst_read_iv(m)) then
+ call cam_pio_get_var(fieldname, fh_ini, arraydims3, arr3d_a, found=readvar)
+ end if
+ call process_inidat('CONSTS', m_cnst=m, fh=fh_ini)
end do
@@ -192,20 +194,19 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
fieldname = 'PHIS'
readvar = .false.
+! if (ideal_phys .or. aqua_planet .or. .not. associated(fh_topo)) then
if (ideal_phys .or. aqua_planet) then
phis_tmp(:,:) = 0._r8
else
- call infld(fieldname, ncid_topo, 'lon', 'lat', 1, plon, 1, plat, &
- phis_tmp, readvar, grid_map='global')
- if(.not. readvar) call endrun('dynamics/eul/inidat.F90')
+ call cam_pio_get_var(fieldname, fh_topo, arraydims2, phis_tmp, found=readvar)
+ if(.not. readvar) call endrun('dynamics/eul/inidat.F90: PHIS not found')
end if
- call process_inidat('PHIS', ncid=ncid_topo)
+ call process_inidat('PHIS', fh=fh_topo)
fieldname = 'PS'
- call infld(fieldname, ncid_ini, 'lon', 'lat', 1, plon, 1, plat, &
- ps_tmp , readvar, grid_map='global')
+ call cam_pio_get_var(fieldname, fh_ini, arraydims2, ps_tmp, found=readvar)
if(.not. readvar) call endrun('PS not found in init file')
@@ -231,8 +232,8 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
if (single_column) then
if ( use_camiop ) then
fieldname = 'CLAT1'
- call infld(fieldname, ncid_ini, 'lon', 'lat', 1, pcols, begchunk, endchunk, &
- clat2d, readvar, grid_map='phys')
+ call infld(fieldname, fh_ini, 'lon', 'lat', 1, pcols, begchunk, endchunk, &
+ clat2d, readvar, gridname='physgrid')
if(.not. readvar) then
call endrun('CLAT not on iop initial file')
else
@@ -241,8 +242,8 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
end if
fieldname = 'CLON1'
- call infld(fieldname, ncid_ini, 'lon', 'lat', 1, pcols, begchunk, endchunk, &
- clon2d, readvar, grid_map='phys')
+ call infld(fieldname, fh_ini, 'lon', 'lat', 1, pcols, begchunk, endchunk, &
+ clon2d, readvar, gridname='physgrid')
if(.not. readvar) then
call endrun('CLON not on iop initial file')
else
@@ -254,10 +255,10 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
! Get latdeg/londeg from initial file for bfb calculations
! needed for dyn_grid to determine bounding area and verticies
!
- ierr = pio_inq_dimid (ncid_ini, 'lon' , londimid)
- ierr = pio_inq_dimlen (ncid_ini, londimid, dimlon)
- ierr = pio_inq_dimid (ncid_ini, 'lat' , latdimid)
- ierr = pio_inq_dimlen (ncid_ini, latdimid, dimlat)
+ ierr = pio_inq_dimid (fh_ini, 'lon' , londimid)
+ ierr = pio_inq_dimlen (fh_ini, londimid, dimlon)
+ ierr = pio_inq_dimid (fh_ini, 'lat' , latdimid)
+ ierr = pio_inq_dimlen (fh_ini, latdimid, dimlat)
strt(:)=1
cnt(1)=dimlon
cnt(2)=dimlat
@@ -265,11 +266,11 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in, datatype)
allocate(latiop(dimlat))
allocate(loniop(dimlon))
allocate(tmp2d(dimlon,dimlat))
- ierr = pio_inq_varid (ncid_ini,'CLAT1', varid)
- ierr = pio_get_var(ncid_ini,varid,strt,cnt,tmp2d)
+ ierr = pio_inq_varid (fh_ini,'CLAT1', varid)
+ ierr = pio_get_var(fh_ini,varid,strt,cnt,tmp2d)
latiop(:)=tmp2d(1,:)
- ierr = pio_inq_varid (ncid_ini,'CLON1', varid)
- ierr = pio_get_var(ncid_ini,varid,strt,cnt,tmp2d)
+ ierr = pio_inq_varid (fh_ini,'CLON1', varid)
+ ierr = pio_get_var(fh_ini,varid,strt,cnt,tmp2d)
loniop(:)=tmp2d(:,1)
deallocate(tmp2d)
else
@@ -315,7 +316,7 @@ end subroutine read_inidat
!*********************************************************************
- subroutine process_inidat(fieldname, m_cnst, ncid)
+ subroutine process_inidat(fieldname, m_cnst, fh)
!
!-----------------------------------------------------------------------
!
@@ -338,9 +339,9 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
use tracers , only: tracers_implements_cnst, tracers_init_cnst
use aoa_tracers , only: aoa_tracers_implements_cnst, aoa_tracers_init_cnst
use clubb_intr , only: clubb_implements_cnst, clubb_init_cnst
- use stratiform, only: stratiform_implements_cnst, stratiform_init_cnst
- use microp_driver, only: microp_driver_implements_cnst, microp_driver_init_cnst
- use phys_control, only: phys_getopts
+ use stratiform , only: stratiform_implements_cnst, stratiform_init_cnst
+ use microp_driver,only: microp_driver_implements_cnst, microp_driver_init_cnst
+ use phys_control, only: phys_getopts
use co2_cycle , only: co2_implements_cnst, co2_init_cnst
use unicon_cam, only: unicon_implements_cnst, unicon_init_cnst
#if ( defined SPMD )
@@ -351,16 +352,17 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
!
! Input arguments
!
- character(len=*), intent(in) :: fieldname ! fields to be processed
- integer , intent(in), optional :: m_cnst ! constituent index
- type(file_desc_t), intent(inout), optional :: ncid ! netcdf file handle
+ character(len=*), intent(in) :: fieldname ! fields to be processed
+ integer, intent(in), optional :: m_cnst ! constituent index
+ type(file_desc_t), intent(inout), optional :: fh ! pio file handle
!
!---------------------------Local workspace-----------------------------
!
integer i,j,k,n,lat,irow ! grid and constituent indices
real(r8) pertval ! perturbation value
integer varid ! netCDF variable id
- integer ret, attlen ! netcdf return values
+ integer ret
+ integer(pio_offset_kind) :: attlen ! netcdf return values
logical phis_hires ! true => PHIS came from hi res topo
character*256 text
character*256 trunits ! tracer untis
@@ -485,58 +487,67 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
allocate ( tmp3d_extend(plon,plev,beglat:endlat) )
-!
-! Check that all tracer units are in mass mixing ratios
-!
- if(readvar) then
- ret = pio_inq_varid(NCID,cnst_name(m_cnst), varid)
- ret = pio_get_att(NCID,varid,'units',trunits)
+ if (readvar) then
+ ! Check that all tracer units are in mass mixing ratios
+ ret = pio_inq_varid(fh, cnst_name(m_cnst), varid)
+ ret = pio_get_att(fh, varid, 'units', trunits)
if (trunits(1:5) .ne. 'KG/KG' .and. trunits(1:5) .ne. 'kg/kg') then
call endrun(' '//trim(subname)//' Error: Units for tracer ' &
//trim(cnst_name(m_cnst))//' must be in KG/KG')
end if
-!
-! Constituents not read from initial file are initialized by the package that implements them.
-!
- else
+
+ else
+ ! Constituents not read from initial file are initialized by the package that implements them.
+
if(m_cnst == 1 .and. moist_physics) then
call endrun(' '//trim(subname)//' Error: Q must be on Initial File')
end if
+ if (masterproc) write(iulog,*) 'Warning: Not reading ',cnst_name(m_cnst), ' from IC file.'
+
arr3d_a(:,:,:) = 0._r8
allocate(gcid(plon))
do j=1,plat
gcid(:) = j
- if(masterproc) write(iulog,*) 'Warning: Not reading ',cnst_name(m_cnst), ' from IC file.'
if (microp_driver_implements_cnst(cnst_name(m_cnst))) then
call microp_driver_init_cnst(cnst_name(m_cnst),arr3d_a(:,:,j) , gcid)
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "microp_driver_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "microp_driver_init_cnst"'
else if (clubb_implements_cnst(cnst_name(m_cnst))) then
call clubb_init_cnst(cnst_name(m_cnst), arr3d_a(:,:,j), gcid)
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "clubb_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "clubb_init_cnst"'
else if (stratiform_implements_cnst(cnst_name(m_cnst))) then
call stratiform_init_cnst(cnst_name(m_cnst), arr3d_a(:,:,j), gcid)
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "stratiform_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "stratiform_init_cnst"'
else if (chem_implements_cnst(cnst_name(m_cnst))) then
call chem_init_cnst(cnst_name(m_cnst), arr3d_a(:,:,j), gcid)
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "chem_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "chem_init_cnst"'
else if (tracers_implements_cnst(cnst_name(m_cnst))) then
call tracers_init_cnst(cnst_name(m_cnst), arr3d_a(:,:,j), gcid)
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "tracers_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "tracers_init_cnst"'
else if (aoa_tracers_implements_cnst(cnst_name(m_cnst))) then
call aoa_tracers_init_cnst(cnst_name(m_cnst), arr3d_a(:,:,j), gcid)
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "aoa_tracers_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "aoa_tracers_init_cnst"'
else if (carma_implements_cnst(cnst_name(m_cnst))) then
call carma_init_cnst(cnst_name(m_cnst), arr3d_a(:,:,j), gcid)
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "carma_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "carma_init_cnst"'
else if (co2_implements_cnst(cnst_name(m_cnst))) then
call co2_init_cnst(cnst_name(m_cnst), arr3d_a(:,:,j), gcid)
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "co2_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "co2_init_cnst"'
else if (unicon_implements_cnst(cnst_name(m_cnst))) then
call unicon_init_cnst(cnst_name(m_cnst), arr3d_a(:,:,j), gcid)
- if (masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' initialized by "unicon_init_cnst"'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' initialized by "unicon_init_cnst"'
else
- if(masterproc) write(iulog,*) ' ', cnst_name(m_cnst), ' set to 0.'
+ if (masterproc .and. j==1) write(iulog,*) ' ', trim(cnst_name(m_cnst)),&
+ ' set to 0.'
end if
end do
@@ -546,7 +557,7 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
!$omp parallel do private(lat)
do lat = 1,plat
call qneg3(trim(subname), lat ,nlon(lat),plon ,plev , &
- m_cnst, m_cnst, qmin(m_cnst) ,arr3d_a(1,1,lat))
+ m_cnst, m_cnst, qmin(m_cnst) ,arr3d_a(1,1,lat), .True.)
end do
!
! if "Q", "CLDLIQ", or "CLDICE", save off for later use
@@ -556,7 +567,7 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
#if ( defined SPMD )
numperlat = plnlv
call compute_gsfactors (numperlat, numrecv, numsend, displs)
- call mpiscatterv (arr3d_a , numsend, displs, mpir8, tmp3d_extend(:,:,beglat:endlat) ,numrecv, mpir8,0,mpicom)
+ call mpiscatterv (arr3d_a , numsend, displs, mpir8, tmp3d_extend ,numrecv, mpir8,0,mpicom)
q3(:,:,m_cnst,:,1) = tmp3d_extend(:,:,beglat:endlat)
#else
q3(:,:plev,m_cnst,:,1) = arr3d_a(:plon,:plev,:plat)
@@ -593,8 +604,8 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
#if ( defined SPMD )
numperlat = plon
call compute_gsfactors (numperlat, numrecv, numsend, displs)
- call mpiscatterv (tmp2d_a ,numsend, displs, mpir8,dpsl(:,beglat:endlat) ,numrecv, mpir8,0,mpicom)
- call mpiscatterv (tmp2d_b ,numsend, displs, mpir8,dpsm(:,beglat:endlat) ,numrecv, mpir8,0,mpicom)
+ call mpiscatterv (tmp2d_a ,numsend, displs, mpir8,dpsl ,numrecv, mpir8,0,mpicom)
+ call mpiscatterv (tmp2d_b ,numsend, displs, mpir8,dpsm ,numrecv, mpir8,0,mpicom)
#else
dpsl(:,:) = tmp2d_a(:,:)
dpsm(:,:) = tmp2d_b(:,:)
@@ -612,14 +623,14 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
! Check for presence of 'from_hires' attribute to decide whether to filter
!
if(readvar) then
- ret = pio_inq_varid (ncid, 'PHIS', varid)
+ ret = pio_inq_varid (fh, 'PHIS', varid)
! Allow pio to return errors in case from_hires doesn't exist
- call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
- ret = pio_inq_attlen (ncid, varid, 'from_hires', attlen)
+ call pio_seterrorhandling(fh, PIO_BCAST_ERROR)
+ ret = pio_inq_attlen (fh, varid, 'from_hires', attlen)
if (ret.eq.PIO_NOERR .and. attlen.gt.256) then
call endrun(' '//trim(subname)//' Error: from_hires attribute length is too long')
end if
- ret = pio_get_att(ncid, varid, 'from_hires', text)
+ ret = pio_get_att(fh, varid, 'from_hires', text)
if (ret.eq.PIO_NOERR .and. text(1:4).eq.'true') then
phis_hires = .true.
@@ -629,7 +640,7 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
if(masterproc) write(iulog,*)trim(subname), ': Will not filter input PHIS: attribute ', &
'from_hires is either false or not present'
end if
- call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
+ call pio_seterrorhandling(fh, PIO_INTERNAL_ERROR)
else
phis_hires = .false.
@@ -647,7 +658,7 @@ subroutine process_inidat(fieldname, m_cnst, ncid)
#if ( defined SPMD )
numperlat = plon
call compute_gsfactors (numperlat, numrecv, numsend, displs)
- call mpiscatterv (phis_tmp ,numsend, displs, mpir8,phis (:,beglat:endlat) ,numrecv, mpir8,0,mpicom)
+ call mpiscatterv (phis_tmp ,numsend, displs, mpir8,phis ,numrecv, mpir8,0,mpicom)
#else
!$omp parallel do private(lat)
do lat = 1,plat
diff --git a/components/cam/src/dynamics/eul/interp_mod.F90 b/components/cam/src/dynamics/eul/interp_mod.F90
index 6d58df43d85..a36f01d7315 100644
--- a/components/cam/src/dynamics/eul/interp_mod.F90
+++ b/components/cam/src/dynamics/eul/interp_mod.F90
@@ -1,9 +1,14 @@
module interp_mod
- use shr_kind_mod, only : r8=>shr_kind_r8
+ use shr_kind_mod, only : r8=>shr_kind_r8
use cam_abortutils, only : endrun
+
implicit none
- public get_interp_lat, get_interp_lon, setup_history_interpolation, write_interpolated
- public var_is_vector_uvar, var_is_vector_vvar, latlon_interpolation, add_interp_attributes
+ private
+ save
+
+ public :: setup_history_interpolation
+ public :: set_interp_hfile
+ public :: write_interpolated
interface write_interpolated
module procedure write_interpolated_scalar
@@ -12,29 +17,27 @@ module interp_mod
integer, parameter :: nlat=0, nlon=0
contains
- subroutine add_interp_attributes(file)
- use pio, only : file_desc_t
- type(file_desc_t) :: file
+ subroutine setup_history_interpolation(interp_ok, mtapes, interp_output, &
+ interp_info)
+ use cam_history_support, only: interp_info_t
- call endrun('This routine is a stub, you shouldnt get here')
-
- end subroutine add_interp_attributes
+ ! Dummy arguments
+ logical, intent(inout) :: interp_ok
+ integer, intent(in) :: mtapes
+ logical, intent(in) :: interp_output(:)
+ type(interp_info_t), intent(inout) :: interp_info(:)
-
- subroutine setup_history_interpolation(mtapes)
- integer, intent(in) :: mtapes
- call endrun('This routine is a stub, you shouldnt get here')
+ interp_ok = .false.
end subroutine setup_history_interpolation
- function latlon_interpolation(t)
- integer, intent(in) :: t
- logical :: latlon_interpolation
-
- latlon_interpolation = .false.
- end function latlon_interpolation
-
+ subroutine set_interp_hfile(hfilenum, interp_info)
+ use cam_history_support, only: interp_info_t
+ ! Dummy arguments
+ integer, intent(in) :: hfilenum
+ type(interp_info_t), intent(inout) :: interp_info(:)
+ end subroutine set_interp_hfile
subroutine write_interpolated_scalar(File, varid, fld, numlev, data_type, decomp_type)
use pio, only : file_desc_t, var_desc_t
@@ -48,9 +51,6 @@ subroutine write_interpolated_scalar(File, varid, fld, numlev, data_type, decomp
end subroutine write_interpolated_scalar
-
-
-
subroutine write_interpolated_vector(File, varidu, varidv, fldu, fldv, numlev, data_type, decomp_type)
use pio, only : file_desc_t, var_desc_t
implicit none
@@ -62,29 +62,4 @@ subroutine write_interpolated_vector(File, varidu, varidv, fldu, fldv, numlev, d
end subroutine write_interpolated_vector
- function var_is_vector_uvar(name)
- character(len=*), intent(in) :: name
- integer ::var_is_vector_uvar
- var_is_vector_uvar=0
- end function var_is_vector_uvar
- function var_is_vector_vvar(name)
- character(len=*), intent(in) :: name
- integer ::var_is_vector_vvar
- var_is_vector_vvar=0
- end function var_is_vector_vvar
-
- function get_interp_lat() result(thislat)
- real(kind=r8) :: thislat(nlat)
- call endrun('This routine is a stub, you shouldnt get here')
-
- return
- end function get_interp_lat
- function get_interp_lon() result(thislon)
- real(kind=r8) :: thislon(nlon)
- call endrun('This routine is a stub, you shouldnt get here')
-
- return
- end function get_interp_lon
-
-
end module interp_mod
diff --git a/components/cam/src/dynamics/eul/restart_dynamics.F90 b/components/cam/src/dynamics/eul/restart_dynamics.F90
index 019413f15df..51d132f57a3 100644
--- a/components/cam/src/dynamics/eul/restart_dynamics.F90
+++ b/components/cam/src/dynamics/eul/restart_dynamics.F90
@@ -3,7 +3,7 @@ module restart_dynamics
use shr_kind_mod, only: r8 => shr_kind_r8
use pio, only : var_desc_t, file_desc_t, pio_double, pio_unlimited, pio_def_var, &
pio_def_dim, io_desc_t, pio_offset_kind, pio_put_var, pio_write_darray, &
- pio_setdebuglevel, pio_setframe, pio_initdecomp, pio_freedecomp, &
+ pio_setdebuglevel,pio_setframe, pio_initdecomp, pio_freedecomp, &
pio_read_darray, pio_inq_varid, pio_get_var
use prognostics, only: u3, v3, t3, q3, &
pdeld, ps, vort, div, &
@@ -168,33 +168,36 @@ end subroutine init_restart_varlist
-subroutine init_restart_dynamics(File, hdimids, dyn_out)
+subroutine init_restart_dynamics(File, dyn_out)
- use dyn_comp, only: dyn_export_t
- use dyn_grid, only: get_horiz_grid_dim_d
- use constituents, only: pcnst
- use hycoef, only: init_restart_hycoef
+ use dyn_comp, only: dyn_export_t
+ use constituents, only: pcnst
+ use hycoef, only: init_restart_hycoef
+ use cam_grid_support, only: cam_grid_write_attr, cam_grid_id
+ use cam_grid_support, only: cam_grid_header_info_t
! Input arguments
type(File_desc_t), intent(inout) :: File
- integer, pointer :: hdimids(:)
type(Dyn_export_t), intent(in) :: dyn_out
+ integer :: hdimids(2)
integer :: vdimids(2)
character(len=namlen) :: name
integer :: alldims(4), alldims2d(3), qdims(5)
- integer :: timelevels_dimid, i, hdim1, hdim2, ierr
+ integer :: timelevels_dimid, i, ierr
type(var_desc_t), pointer :: vdesc
+ integer :: grid_id
integer :: ndims, timelevels
+ type(cam_grid_header_info_t) :: info
call init_restart_hycoef(File, vdimids)
- call get_horiz_grid_dim_d(hdim1, hdim2)
- allocate(hdimids(2))
- ierr = PIO_Def_Dim(File, 'lon',hdim1, hdimids(1))
-
- ierr = PIO_Def_Dim(File, 'lat',hdim2, hdimids(2))
+ ! Grid attributes
+ grid_id = cam_grid_id('gauss_grid')
+ call cam_grid_write_attr(File, grid_id, info)
+ hdimids(1) = info%get_hdimid(1)
+ hdimids(2) = info%get_hdimid(2)
ierr = PIO_Def_Dim(File,'timelevels',PIO_UNLIMITED,timelevels_dimid)
@@ -227,7 +230,7 @@ subroutine init_restart_dynamics(File, hdimids, dyn_out)
do i=1,restartvarcnt
- call get_restart_var(File, i, name, timelevels, ndims, vdesc)
+ call get_restart_var(i, name, timelevels, ndims, vdesc)
if(timelevels>1) then
if(ndims==3) then
ierr = PIO_Def_Var(File, name, pio_double, alldims2d, vdesc)
@@ -265,6 +268,8 @@ subroutine write_restart_dynamics (File, dyn_out)
use constituents, only: pcnst
use eul_control_mod, only: fixmas, tmass0
use hycoef, only: write_restart_hycoef
+ use cam_grid_support, only: cam_grid_write_var
+ use dyn_grid, only: dyn_decomp
!
@@ -288,6 +293,9 @@ subroutine write_restart_dynamics (File, dyn_out)
character(len=namlen) :: name
!
+ ! Write grid vars
+ call cam_grid_write_var(File, dyn_decomp)
+
call write_restart_hycoef(File)
call get_curr_time(ndcur, nscur)
@@ -319,7 +327,7 @@ subroutine write_restart_dynamics (File, dyn_out)
ierr = pio_put_var(File,timedesc%varid, (/int(t)/), time)
end do
do i=1,restartvarcnt
- call get_restart_var(File, i, name, timelevels, ndims, vdesc)
+ call get_restart_var(i, name, timelevels, ndims, vdesc)
if(timelevels==1) then
if(ndims==2) then
call pio_write_darray(File, vdesc, iodesc2d, transfer(restartvars(i)%v2d(:,:), mold), ierr)
@@ -334,7 +342,7 @@ subroutine write_restart_dynamics (File, dyn_out)
if(t==2) ct=n3m1
if(t==3) ct=n3
- call pio_setframe(vdesc, t)
+ call pio_setframe(File, vdesc, t)
if(ndims==3) then
call pio_write_darray(File, vdesc, iodesc2d, transfer(restartvars(i)%v3d(:,:,ct), mold), ierr)
else if(ndims==4) then
@@ -354,8 +362,7 @@ subroutine write_restart_dynamics (File, dyn_out)
return
end subroutine write_restart_dynamics
- subroutine get_restart_var(File, i,name, timelevels, ndims, vdesc)
- type(file_desc_t), intent(in) :: File ! PIO file handle
+ subroutine get_restart_var(i,name, timelevels, ndims, vdesc)
integer, intent(in) :: i
character(len=namlen), intent(out) :: name
integer, intent(out) :: ndims, timelevels
@@ -368,7 +375,6 @@ subroutine get_restart_var(File, i,name, timelevels, ndims, vdesc)
allocate(restartvars(i)%vdesc)
end if
vdesc => restartvars(i)%vdesc
- call pio_setframe(File, vdesc, int(-1,pio_offset_kind))
end subroutine get_restart_var
@@ -462,7 +468,7 @@ subroutine read_restart_dynamics (File, dyn_in, dyn_out, NLFileName)
call init_iop_fields()
#endif
do i=1,restartvarcnt
- call get_restart_var(File, i, name, timelevels, ndims, vdesc)
+ call get_restart_var(i, name, timelevels, ndims, vdesc)
ierr = PIO_Inq_varid(File, name, vdesc)
@@ -483,7 +489,7 @@ subroutine read_restart_dynamics (File, dyn_in, dyn_out, NLFileName)
if(t==1) ct=n3m2
if(t==2) ct=n3m1
if(t==3) ct=n3
- call pio_setframe(vdesc, t)
+ call pio_setframe(File, vdesc, t)
if(ndims==3) then
call pio_read_darray(File, vdesc, iodesc2d, tmp(1:s2d), ierr)
restartvars(i)%v3d(:,:,ct) = reshape(tmp(1:s2d), dims2d)
diff --git a/components/cam/src/dynamics/eul/scanslt.F90 b/components/cam/src/dynamics/eul/scanslt.F90
index e0391766b87..209c64c238a 100644
--- a/components/cam/src/dynamics/eul/scanslt.F90
+++ b/components/cam/src/dynamics/eul/scanslt.F90
@@ -9,11 +9,11 @@ module scanslt
! $Id$
!
!-----------------------------------------------------------------------
- use shr_kind_mod, only: r8 => shr_kind_r8
- use pmgrid, only: plon, plat, plev, beglat, endlat, plevp
- use constituents, only: pcnst
+ use shr_kind_mod, only: r8 => shr_kind_r8
+ use pmgrid, only: plon, plat, plev, beglat, endlat, plevp
+ use constituents, only: pcnst
use cam_abortutils, only: endrun
- use scamMod, only: single_column
+ use scamMod, only: single_column
use perf_mod
!-----------------------------------------------------------------------
implicit none
@@ -131,7 +131,7 @@ end subroutine scanslt_alloc
!
!-----------------------------------------------------------------------
!
-subroutine scanslt_initial( adv_state, etamid, gravit_in, gw, detam, cwava )
+subroutine scanslt_initial( adv_state, etamid, gravit_in, detam, cwava )
!-----------------------------------------------------------------------
!
! Purpose:
@@ -146,17 +146,16 @@ subroutine scanslt_initial( adv_state, etamid, gravit_in, gw, detam, cwava )
use prognostics, only: ps, n3
use rgrid, only: nlon
use time_manager, only: is_first_step
- use hycoef, only: hyai, hybi, ps0
+ use hycoef, only: hyam, hybm, hyai, hybi, ps0
use eul_control_mod, only : pdela
!
! Input arguments
!
- real(r8), intent(inout) :: etamid(plev) ! vertical coords at midpoints
+ real(r8), intent(out) :: etamid(plev) ! vertical coords at midpoints
real(r8), intent(in) :: gravit_in ! Gravitational constant
!
! Output arguments
!
- real(r8), intent(out) :: gw(plat) ! Gaussian weights
real(r8), intent(out) :: detam(plev) ! intervals between vert full levs.
real(r8), intent(out) :: cwava(plat) ! weight applied to global integrals
type(advection_state), intent(out) :: adv_state ! Advection state data
@@ -169,16 +168,17 @@ subroutine scanslt_initial( adv_state, etamid, gravit_in, gw, detam, cwava )
real(r8) :: pmid(plon,plev) ! pressure at model levels
real(r8) :: pint(plon,plevp) ! pressure at interfaces
real(r8) :: pdel(plon,plev) ! pressure difference between
+ real(r8) :: gw(plat) ! Gaussian weights needed for SCAM grdini call
!
! Allocate memory for scanslt variables
!
call adv_state_alloc( adv_state )
-!
-! Eta at interfaces
-!
- do k=1,plevp
+
+ do k = 1, plev
+ etamid(k) = hyam(k) + hybm(k)
etaint(k) = hyai(k) + hybi(k)
end do
+ etaint(plevp) = hyai(plevp) + hybi(plevp)
!
! For SCAM compute pressure levels to use for eta interface
!
@@ -636,7 +636,8 @@ subroutine grdini(pmap ,etamid ,etaint ,gravit ,dlam , &
! Reviewed: D. Williamson, P. Rasch, March 1996
!
!-----------------------------------------------------------------------
- use rgrid, only: nlon
+ use rgrid, only: nlon
+ use vrtmap_mod, only: vrtmap
!------------------------------Parameters-------------------------------
!
! Input arguments
@@ -716,48 +717,48 @@ subroutine grdini(pmap ,etamid ,etaint ,gravit ,dlam , &
!-----------------------------------------------------------------------
if (single_column) then
- dlam(:)=0._r8
- lam(:,:)=0._r8
- phi(:)=0._r8
- dphi(:)=0._r8
- sinlam(:,:)=0._r8
- coslam(:,:)=0._r8
- detai(:)=0._r8
- kdpmpf(:)=0._r8
- kdpmph(:)=0._r8
- gw(:)=1._r8
- call basdz(plev ,etamid ,lbasdz )
- call basdz(plevp ,etaint ,lbassd )
+ dlam(:)=0._r8
+ lam(:,:)=0._r8
+ phi(:)=0._r8
+ dphi(:)=0._r8
+ sinlam(:,:)=0._r8
+ coslam(:,:)=0._r8
+ detai(:)=0._r8
+ kdpmpf(:)=0._r8
+ kdpmph(:)=0._r8
+ gw(:)=1._r8
+ call basdz(plev ,etamid ,lbasdz )
+ call basdz(plevp ,etaint ,lbassd )
else
-!
-! Initialize extended horizontal grid coordinates.
-!
- call grdxy(dlam ,lam ,phi ,gw ,sinlam , &
- coslam )
-!
-! Basis functions for computing Lagrangian cubic derivatives
-! on unequally spaced latitude and vertical grids.
-!
- call basdy(phi ,lbasdy )
-
- call basdz(plev ,etamid ,lbasdz )
- call basdz(plevp ,etaint ,lbassd )
-
-
-!
-! Basis functions for computing weights for Lagrangian cubic
-! interpolation on unequally spaced latitude grids.
-!
- call basiy(phi ,lbasiy )
-!
-! Compute interval lengths in latitudinal grid
-!
- do j = 1,platd-1
- dphi(j) = phi(j+1) - phi(j)
- end do
+ !
+ ! Initialize extended horizontal grid coordinates.
+ !
+ call grdxy(dlam ,lam ,phi ,gw ,sinlam , &
+ coslam )
+ !
+ ! Basis functions for computing Lagrangian cubic derivatives
+ ! on unequally spaced latitude and vertical grids.
+ !
+ call basdy(phi ,lbasdy )
+
+ call basdz(plev ,etamid ,lbasdz )
+ call basdz(plevp ,etaint ,lbassd )
+
+
+ !
+ ! Basis functions for computing weights for Lagrangian cubic
+ ! interpolation on unequally spaced latitude grids.
+ !
+ call basiy(phi ,lbasiy )
+ !
+ ! Compute interval lengths in latitudinal grid
+ !
+ do j = 1,platd-1
+ dphi(j) = phi(j+1) - phi(j)
+ end do
-endif
+ endif
!
! Compute interval lengths in vertical grids.
!
diff --git a/components/cam/src/dynamics/eul/stepon.F90 b/components/cam/src/dynamics/eul/stepon.F90
index 1bc00788960..e21948ba35d 100644
--- a/components/cam/src/dynamics/eul/stepon.F90
+++ b/components/cam/src/dynamics/eul/stepon.F90
@@ -15,16 +15,15 @@ module stepon
use camsrfexch, only: cam_out_t
use ppgrid, only: begchunk, endchunk
use physics_types, only: physics_state, physics_tend
- use time_manager, only: is_first_step
+ use time_manager, only: is_first_step, get_step_size
use iop, only: setiopupdate, readiopdata
use scamMod, only: use_iop,doiopupdate,use_pert_frc,wfld,wfldh,single_column
use perf_mod
+
implicit none
+ private
+ save
- private ! By default make all data and methods private to this module
-!
-! Public methods
-!
public stepon_init ! Initialization
public stepon_run1 ! Run method phase 1
public stepon_run2 ! Run method phase 2
@@ -33,7 +32,6 @@ module stepon
!
! Private module data
!
- save
type(physics_state), pointer :: phys_state(:) ! Physics state data
type(physics_tend ), pointer :: phys_tend(:) ! Physics tendency data
@@ -53,18 +51,15 @@ module stepon
real(r8) :: pdel(plon,plev) ! Pressure depth of layer
real(r8) :: pint(plon,plevp) ! Pressure at interfaces
real(r8) :: pmid(plon,plev) ! Pressure at midpoint
- real(r8) :: dtime ! timestep size
type(advection_state) :: adv_state ! Advection state data
+ real(r8) :: etamid(plev) ! vertical coords at midpoints or pmid if single_column
+
!=======================================================================
contains
!=======================================================================
-!
-!=======================================================================
-!
-
-subroutine stepon_init( gw, etamid, dyn_in, dyn_out )
+subroutine stepon_init(dyn_in, dyn_out)
!-----------------------------------------------------------------------
!
! Purpose: Initialization, primarily of dynamics.
@@ -76,8 +71,6 @@ subroutine stepon_init( gw, etamid, dyn_in, dyn_out )
use constituents, only: pcnst
use physconst, only: gravit
use rgrid, only: nlon
- use hycoef, only: hyam, hybm
- use time_manager, only: get_step_size
use eul_control_mod,only: eul_nsplit
#if ( defined BFB_CAM_SCAM_IOP )
use iop, only:init_iop_fields
@@ -85,8 +78,6 @@ subroutine stepon_init( gw, etamid, dyn_in, dyn_out )
!-----------------------------------------------------------------------
! Arguments
!
- real(r8), intent(out) :: gw(plat) ! Gaussian weights
- real(r8), intent(out) :: etamid(plev) ! vertical coords at midpoints
type(dyn_import_t) :: dyn_in ! included for compatibility
type(dyn_export_t) :: dyn_out ! included for compatibility
!-----------------------------------------------------------------------
@@ -97,16 +88,7 @@ subroutine stepon_init( gw, etamid, dyn_in, dyn_out )
call t_startf ('stepon_startup')
- dtime = get_step_size()
- !
- ! Define eta coordinates: Used for calculation etadot vertical velocity
- ! for slt.
- !
- do k=1,plev
- etamid(k) = hyam(k) + hybm(k)
- end do
-
- call scanslt_initial( adv_state, etamid, gravit, gw, detam, cwava )
+ call scanslt_initial(adv_state, etamid, gravit, detam, cwava)
!
! Initial guess for trajectory midpoints in spherical coords.
! nstep = 0: use arrival points as initial guess for trajectory midpoints.
@@ -195,25 +177,20 @@ subroutine stepon_run1( ztodt, phys_state, phys_tend , pbuf2d, dyn_in, dyn_out)
type(dyn_import_t) :: dyn_in ! included for compatibility
type(dyn_export_t) :: dyn_out ! included for compatibility
+ real(r8) :: dtime ! timestep size
+ !-----------------------------------------------------------------------
+
+ dtime = get_step_size()
+
ztodt = 2.0_r8*dtime
- !
+
! If initial time step adjust dt
- !
if (is_first_step()) ztodt = dtime
! subcycling case, physics dt is always dtime
if (eul_nsplit>1) ztodt = dtime
- !
- ! adjust hydrostatic matrices if the time step has changed. This only
- ! happens on transition from time 0 to time 1.
- ! (this is now initialized in dynpgk)
- !if (get_nstep() == 1) then
- ! call settau(dtime)
- !end if
- !
! Dump state variables to IC file
- !
call t_startf ('diag_dynvar_ic')
call diag_dynvar_ic (phis, ps(:,beglat:endlat,n3m1), t3(:,:,beglat:endlat,n3m1), u3(:,:,beglat:endlat,n3m1), &
v3(:,:,beglat:endlat,n3m1), q3(:,:,:,beglat:endlat,n3m1) )
@@ -250,7 +227,7 @@ subroutine stepon_run2( phys_state, phys_tend, dyn_in, dyn_out )
call t_startf ('p_d_coupling')
call p_d_coupling (phys_state, phys_tend, t2, fu, fv, flx_net, &
- qminus(:,:,:,beglat:endlat) )
+ qminus )
call t_stopf ('p_d_coupling')
end subroutine stepon_run2
@@ -258,7 +235,7 @@ end subroutine stepon_run2
!=======================================================================
!
-subroutine stepon_run3( ztodt, etamid, cam_out, phys_state, dyn_in, dyn_out )
+subroutine stepon_run3( ztodt, cam_out, phys_state, dyn_in, dyn_out )
!-----------------------------------------------------------------------
!
! Purpose: Final phase of dynamics run method. Run the actual dynamics.
@@ -268,7 +245,6 @@ subroutine stepon_run3( ztodt, etamid, cam_out, phys_state, dyn_in, dyn_out )
use eul_control_mod,only: eul_nsplit
real(r8), intent(in) :: ztodt ! twice time step unless nstep=0
type(cam_out_t), intent(inout) :: cam_out(begchunk:endchunk)
- real(r8), intent(in) :: etamid(plev) ! vertical coords at midpoints
type(physics_state), intent(in):: phys_state(begchunk:endchunk)
type(dyn_import_t) :: dyn_in ! included for compatibility
type(dyn_export_t) :: dyn_out ! included for compatibility
diff --git a/components/cam/src/dynamics/eul/tfilt_massfix.F90 b/components/cam/src/dynamics/eul/tfilt_massfix.F90
index 952a43d7870..76d20f080fb 100644
--- a/components/cam/src/dynamics/eul/tfilt_massfix.F90
+++ b/components/cam/src/dynamics/eul/tfilt_massfix.F90
@@ -316,7 +316,7 @@ subroutine tfilt_massfixrun (ztodt, lat, u3m1, u3, &
! Check for and correct invalid constituents
!
call qneg3 ('TFILT_MASSFIX',lat ,nlon ,plon ,plev , &
- 1, pcnst, qmin ,q3(1,1,1))
+ 1, pcnst, qmin ,q3(1,1,1),.True.)
!
! Send slt tendencies to the history tape
!
diff --git a/components/cam/src/utils/spmd_utils.F90 b/components/cam/src/utils/spmd_utils.F90
index 59249244fa0..d0b6df2308a 100644
--- a/components/cam/src/utils/spmd_utils.F90
+++ b/components/cam/src/utils/spmd_utils.F90
@@ -52,7 +52,8 @@ module spmd_utils
! Forward from mpishorthand.F with the idea of phasing out use of and removing that file
!
#ifndef SPMD
- integer :: mpi_status_ignore ! Needs to be defined in mpi-serial
+! MPI_STATUS_IGNORE(MPI_STATUS_SIZE) defined in mpi-serial/mpif.h. commented out by wlin for SCM
+! integer :: mpi_status_ignore ! Needs to be defined in mpi-serial
integer :: mpir8
#endif
!
diff --git a/components/cam/src/utils/vrtmap.F90 b/components/cam/src/utils/vrtmap.F90
new file mode 100644
index 00000000000..5097287e9c5
--- /dev/null
+++ b/components/cam/src/utils/vrtmap.F90
@@ -0,0 +1,77 @@
+module vrtmap_mod
+
+implicit none
+
+private
+
+public :: vrtmap
+
+contains
+
+subroutine vrtmap (pkdim ,pmap ,sigln ,dsigln ,kdpmap )
+!-----------------------------------------------------------------------
+!
+! Purpose: Map indices of an artificial evenly spaced (in log) vertical grid to
+! the indices of the log of the model vertical grid. The resultant
+! array of mapped indices will be used by "kdpfnd" to find the vertical
+! location of any departure point relative to the model grid.
+!
+! Method:
+!
+! Author: Jerry Olson
+!
+!-----------------------------------------------------------------------
+ use shr_kind_mod, only: r8 => shr_kind_r8
+ use cam_abortutils, only: endrun
+ use cam_logfile, only: iulog
+!-----------------------------------------------------------------------
+!
+! Arguments
+!
+ integer, intent(in) :: pkdim ! dimension of "sigln" and "dsigln"
+ integer, intent(in) :: pmap ! dimension of "kdpmap"
+
+ real(r8), intent(in) :: sigln (pkdim) ! model levels (log(eta))
+ real(r8), intent(in) :: dsigln(pkdim) ! intervals between model levels (log)
+
+ integer, intent(out) :: kdpmap(pmap) ! array of mapped indices
+!
+!---------------------------Local variables-----------------------------
+!
+ integer imin ! |
+ integer k ! |-- indices
+ integer kk ! |
+ integer newmap ! estimated value of "pmap"
+
+ real(r8) del ! artificial grid interval
+ real(r8) dp ! artificial departure point
+ real(r8) eps ! epsilon factor
+!
+!-----------------------------------------------------------------------
+!
+ eps = 1.e-05_r8
+ del = ( sigln(pkdim) - sigln(1) )/real(pmap,r8)
+ imin = minloc( dsigln(:pkdim-1), dim=1 )
+ if (del + eps >= dsigln(imin)) then
+ newmap = ( sigln(pkdim) - sigln(1) )/dsigln(imin) + 1
+ write(iulog,9000) pmap,newmap
+ call endrun()
+ end if
+
+ kdpmap(1) = 1
+ do kk = 2,pmap
+ dp = sigln(1) + real(kk-1,r8)*del
+ do k = 1,pkdim-1
+ if(dp > sigln(k) + eps) then
+ kdpmap(kk) = k
+ end if
+ end do
+ end do
+
+ return
+9000 format(' VRTMAP: Not enough artificial grid intervals.'/ &
+ ' Currently, "pmap" is set to ',i20/ &
+ ' Reset parameter "pmap" to at least ',i20)
+end subroutine vrtmap
+
+end module vrtmap_mod