<p><b>duda</b> 2012-06-21 14:01:09 -0600 (Thu, 21 Jun 2012)</p><p>BRANCH COMMIT<br>
<br>
- update the MPAS_TO_CAM() interface to provide 3d pressure and height at layer<br>
centers and interfaces to CAM, since we can no longer compute 3d pressure from<br>
surface pressure and A and B arrays.<br>
<br>
- add basic implementations of interface routines to support a basic cold start<br>
run of MPAS-A from within CAM.<br>
<br>
<br>
M src/driver/mpas_cam_interface.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/cam_mpas_nh/src/driver/mpas_cam_interface.F
===================================================================
--- branches/cam_mpas_nh/src/driver/mpas_cam_interface.F 2012-06-21 18:27:14 UTC (rev 1999)
+++ branches/cam_mpas_nh/src/driver/mpas_cam_interface.F 2012-06-21 20:01:09 UTC (rev 2000)
@@ -6,6 +6,7 @@
use mpas_io_output
use mpas_dmpar
use mpas_constants
+ use mpas_core
type (dm_info), pointer :: dminfo
@@ -14,6 +15,7 @@
integer :: itimestep, n_subcycle_steps
type (io_output_object), save :: output_obj
logical :: restart_run
+ integer :: output_frame
contains
@@ -43,13 +45,9 @@
write(0,*) 'Called MPAS_INIT1'
-#if 0
- allocate(dminfo)
- call dmpar_init(dminfo, mpi_comm)
+ call mpas_framework_init(dminfo, domain, mpi_comm)
- call read_namelist(dminfo)
-
- itimestep = 0
+ itimestep = 1
restart_run = .false. ! If a restart run, this will be set to true in cam_restart_to_mpas
!
@@ -60,12 +58,11 @@
n_subcycle_steps = ceiling(dt_physics / config_dt)
- call allocate_domain(domain, dminfo)
- ! NOTE: Really, we only want to read in the grid information, since any state
- ! variables will be overwritten by cam_inidat_to_mpas() or cam_restart_to_mpas();
- ! however, with a suitable registry, this shouldn't be a problem
- call input_state_for_domain(domain)
+ !
+ ! Read initial conditions
+ !
+ call mpas_input_state_for_domain(domain)
latCell => domain % blocklist % mesh % latCell % array
lonCell => domain % blocklist % mesh % lonCell % array
@@ -89,7 +86,6 @@
call r3_normalize(north(1,iCell), north(2,iCell), north(3,iCell))
end do
-#endif
end subroutine mpas_init1
@@ -142,17 +138,18 @@
write(0,*) 'Called CAM_INIDAT_TO_MPAS'
+ p0 = Pref !MGD this also appears below...
#if 0
block => domain % blocklist
- surface_pressure => block % time_levs(1) % state % surface_pressure % array
- t_pot => block % time_levs(1) % state % theta % array
- u_normal => block % time_levs(1) % state % u % array
- pressure => block % time_levs(1) % state % pressure % array
- geopotential => block % time_levs(1) % state % geopotential % array
- scalars => block % time_levs(1) % state % scalars % array
- alpha => block % time_levs(1) % state % alpha % array
- h => block % time_levs(1) % state % h % array
+ surface_pressure => block % state % time_levs(1) % state % surface_pressure % array
+ t_pot => block % state % time_levs(1) % state % theta % array
+ u_normal => block % state % time_levs(1) % state % u % array
+ pressure => block % state % time_levs(1) % state % pressure % array
+ geopotential => block % state % time_levs(1) % state % geopotential % array
+ scalars => block % state % time_levs(1) % state % scalars % array
+ alpha => block % state % time_levs(1) % state % alpha % array
+ h => block % state % time_levs(1) % state % h % array
rdnu => block % mesh % rdnu % array
rdnw => block % mesh % rdnw % array
fnm => block % mesh % fnm % array
@@ -173,8 +170,8 @@
write(0,*) 'Error: mismatch between Plev and nVertLevels: ', Plev, block % mesh % nVertLevels
return
end if
- if (Pcnst /= num_scalars) then
- write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, num_scalars
+ if (Pcnst /= block % state % time_levs(1) % state % num_scalars) then
+ write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, block % state % time_levs(1) % state % num_scalars
return
end if
@@ -319,29 +316,29 @@
!
! Do a ghost-cell update
!
- call dmpar_exch_halo_field1dReal(domain % dminfo, block % time_levs(1) % state % surface_pressure % array(:), &
- block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % theta % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % pressure % array(:,:), &
- block % mesh % nVertLevels+1, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % geopotential % array(:,:), &
- block % mesh % nVertLevels+1, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(1) % state % scalars % array(:,:,:), &
- num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call mpas_dmpar_exch_halo_field(domain % dminfo, block % state % time_levs(1) % state % surface_pressure % array(:), &
+ block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call mpas_dmpar_exch_halo_field(domain % dminfo, block % state % time_levs(1) % state % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call mpas_dmpar_exch_halo_field(domain % dminfo, block % state % time_levs(1) % state % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call mpas_dmpar_exch_halo_field(domain % dminfo, block % state % time_levs(1) % state % pressure % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call mpas_dmpar_exch_halo_field(domain % dminfo, block % state % time_levs(1) % state % geopotential % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call mpas_dmpar_exch_halo_field(domain % dminfo, block % state % time_levs(1) % state % scalars % array(:,:,:), &
+ block % state % time_levs(1) % state % num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
do i=2,nTimeLevs
- call copy_state(block % time_levs(1) % state, block % time_levs(i) % state)
+ call copy_state(block % state % time_levs(1) % state, block % state % time_levs(i) % state)
end do
- block % time_levs(1) % state % xtime % scalar = 0.0
+ block % state % time_levs(1) % state % xtime % scalar = 0.0
#endif
end subroutine cam_inidat_to_mpas
@@ -400,37 +397,37 @@
write(0,*) 'Error: mismatch between Plev and nVertLevels: ', Plev, block % mesh % nVertLevels
return
end if
- if (Pcnst /= num_scalars) then
- write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, num_scalars
+ if (Pcnst /= block % state % time_levs(1) % state % num_scalars) then
+ write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, block % state % time_levs(1) % state % num_scalars
return
end if
!
- ! Fill in output arrays from block % time_levs(1) % state arrays
+ ! Fill in output arrays from block % state % time_levs(1) % state arrays
! NOTE: We must copy from time_levs(1) rather than time_levs(2), since
! mpas_dyn_run swaps the new time level into time_levs(1) after
! taking a timestep.
!
do iCell=1,block % mesh % nCellsSolve
- Psd(iCell) = block % time_levs(1) % state % surface_pressure % array(iCell)
- Phis(iCell) = block % time_levs(1) % state % geopotential % array(1,iCell)
+ Psd(iCell) = block % state % time_levs(1) % state % surface_pressure % array(iCell)
+ Phis(iCell) = block % state % time_levs(1) % state % geopotential % array(1,iCell)
do k=1,block % mesh % nVertLevels
- Theta(iCell,k) = block % time_levs(1) % state % theta % array(k,iCell)
- W(iCell,k) = block % time_levs(1) % state % ww % array(k,iCell)
+ Theta(iCell,k) = block % state % time_levs(1) % state % theta % array(k,iCell)
+ W(iCell,k) = block % state % time_levs(1) % state % ww % array(k,iCell)
!
! NOTE: Eventually, we need to ensure that the moisture variables we return to CAM
! are what CAM expects, rather than simply what we have in our scalars array
!
- do iScalar=1,num_scalars
- Tracer(iCell,k,iScalar) = block % time_levs(1) % state % scalars % array(iScalar,k,iCell)
+ do iScalar=1,block % state % time_levs(1) % state % num_scalars
+ Tracer(iCell,k,iScalar) = block % state % time_levs(1) % state % scalars % array(iScalar,k,iCell)
end do
do i=1,nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
- U(iCell,i,k) = block % time_levs(1) % state % u % array(k,iEdge)
+ U(iCell,i,k) = block % state % time_levs(1) % state % u % array(k,iEdge)
end do
end do
- W(iCell,k) = block % time_levs(1) % state % ww % array(k,iCell)
+ W(iCell,k) = block % state % time_levs(1) % state % ww % array(k,iCell)
end do
#endif
@@ -496,14 +493,14 @@
write(0,*) 'Error: mismatch between Plev and nVertLevels: ', Plev, block % mesh % nVertLevels
return
end if
- if (Pcnst /= num_scalars) then
- write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, num_scalars
+ if (Pcnst /= block % state % time_levs(1) % state % num_scalars) then
+ write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, block % state % time_levs(1) % state % num_scalars
return
end if
!
- ! Fill in MPAS fields in block % time_levs(1) % state arrays from arguments
+ ! Fill in MPAS fields in block % state % time_levs(1) % state arrays from arguments
! Initialize most variables using cam_inidat_to_mpas
!
call cam_inidat_to_mpas(Numcols, MaxEdges, Plev, Pcnst, Pref, Ptop, &
@@ -511,7 +508,7 @@
do iCell=1,Numcols
do k=1,Plev+1
- block % time_levs(1) % state % ww % array(k,iCell) = W(iCell,k)
+ block % state % time_levs(1) % state % ww % array(k,iCell) = W(iCell,k)
end do
end do
@@ -519,12 +516,12 @@
!
! Do a ghost-cell update
!
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % ww % array(:,:), &
+ call mpas_dmpar_exch_halo_field(domain % dminfo, block % state % time_levs(1) % state % ww % array(:,:), &
block % mesh % nVertLevels+1, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
do i=2,nTimeLevs
- block % time_levs(i) % state % ww % array(:,:) = block % time_levs(1) % state % ww % array(:,:)
+ block % state % time_levs(i) % state % ww % array(:,:) = block % state % time_levs(1) % state % ww % array(:,:)
end do
#endif
@@ -536,7 +533,7 @@
! SUBROUTINE MPAS_INIT2
!
! Calls the core-specific (most likely the hydrostatic atmosphere core)
- ! initialization routine mpas_init() after initial fields have been provided
+ ! initialization routine mpas_core_init() after initial fields have been provided
! through a call to either cam_inidat_to_mpas() or cam_restart_to_mpas().
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -546,6 +543,9 @@
integer :: i, iCell, k, iScalar
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
type (block_type), pointer :: block
+
+ integer :: index_qv
+ character (len=StrKIND) :: initialTimeStamp
write(0,*) 'Called MPAS_INIT2'
@@ -556,42 +556,50 @@
!
if (.not. restart_run) then
block => domain % blocklist
+
+ index_qv = block % state % time_levs(1) % state % index_qv
+
do while (associated(block))
- scalars => block % time_levs(1) % state % scalars % array
+ scalars => block % state % time_levs(1) % state % scalars % array
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
- do iScalar=1,num_scalars
+ do iScalar=1,block % state % time_levs(1) % state % num_scalars
scalars(iScalar,k,iCell) = scalars(iScalar,k,iCell) * (1.0 + scalars(index_qv,k,iCell))
end do
end do
end do
do i=2,nTimeLevs
- block % time_levs(i) % state % scalars % array(:,:,:) = block % time_levs(1) % state % scalars % array(:,:,:)
+ block % state % time_levs(i) % state % scalars % array(:,:,:) = block % state % time_levs(1) % state % scalars % array(:,:,:)
end do
block => block % next
end do
end if
+#endif
! Compute diagnostic fields needed in solve loop, and initialize
! simulation time to 0 for all blocks
block => domain % blocklist
do while (associated(block))
- call mpas_init(block, block % mesh, dt_dynamics)
+ call mpas_core_init(domain, initialTimeStamp)
block => block % next
end do
output_frame = 1
- call output_state_init(output_obj, domain, "OUTPUT")
- call output_state_for_domain(output_obj, domain, output_frame)
-#endif
-
+ if(config_frames_per_outfile > 0) then
+ call mpas_output_state_init(output_obj, domain, "OUTPUT", trim(initialTimeStamp))
+ else
+ call mpas_output_state_init(output_obj, domain, "OUTPUT")
+ end if
+ call atm_write_output_frame(output_obj, output_frame, domain)
+
end subroutine mpas_init2
- subroutine mpas_to_cam(Numcols, Plev, Pcnst, Psd, Phis, T, Ux, Uy, Omega, Tracer)
+ subroutine mpas_to_cam(Numcols, Plev, Pcnst, Psd, Phis, Pint, Pmid, Zint, Zmid, &
+ T, Ux, Uy, Omega, Tracer)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SUBROUTINE MPAS_TO_CAM
!
@@ -601,6 +609,10 @@
!
! Output: Psd - dry surface pressure (Pa)
! Phis - surface geopotential (m^2/s^2)
+ ! Pint - dry pressure at vertical layer interfaces
+ ! Pmid - dry pressure at vertical layer mid-points
+ ! Zint - geopotential height at vertical layer interfaces
+ ! Zmid - geopotential height at vertical layer mid-points
! T - temperature (K)
! Ux - longitudinal velocity at cell centers (m/s)
! Uy - latitudinal velocity at cell centers (m/s)
@@ -616,6 +628,10 @@
integer, intent(in) :: Pcnst
real (kind=RKIND), dimension(Numcols), intent(out) :: Psd
real (kind=RKIND), dimension(Numcols), intent(out) :: Phis
+ real (kind=RKIND), dimension(Numcols,Plev+1), intent(out) :: Pint
+ real (kind=RKIND), dimension(Numcols,Plev), intent(out) :: Pmid
+ real (kind=RKIND), dimension(Numcols,Plev+1), intent(out) :: Zint
+ real (kind=RKIND), dimension(Numcols,Plev), intent(out) :: Zmid
real (kind=RKIND), dimension(Numcols,Plev), intent(out) :: T
real (kind=RKIND), dimension(Numcols,Plev), intent(out) :: Ux
real (kind=RKIND), dimension(Numcols,Plev), intent(out) :: Uy
@@ -623,26 +639,18 @@
real (kind=RKIND), dimension(Numcols,Plev,Pcnst), intent(out) :: Tracer
integer :: iCell, k, iScalar
+ integer :: index_qv
real (kind=RKIND), dimension(:), pointer :: latCell, lonCell
- real (kind=RKIND), dimension(:,:), pointer :: theta, pressure, ww, uReconstZonal, uReconstMeridional, east, north
+ real (kind=RKIND), dimension(:,:), pointer :: theta, pressure_base, pressure_p, ww, uReconstZonal, uReconstMeridional, east, north
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
type (block_type), pointer :: block
write(0,*) 'Called MPAS_TO_CAM'
-#if 0
+
block => domain % blocklist
- latCell => block % mesh % latCell % array
- lonCell => block % mesh % lonCell % array
- east => block % mesh % east % array
- north => block % mesh % north % array
- theta => block % time_levs(1) % state % theta % array
- pressure => block % time_levs(1) % state % pressure % array
- ww => block % time_levs(1) % state % ww % array
- uReconstZonal => block % time_levs(1) % state % uReconstructZonal % array
- uReconstMeridional => block % time_levs(1) % state % uReconstructMeridional % array
- scalars => block % time_levs(1) % state % scalars % array
+ index_qv = block % state % time_levs(1) % state % index_qv
!
@@ -656,29 +664,50 @@
write(0,*) 'Error: mismatch between Plev and nVertLevels: ', Plev, block % mesh % nVertLevels
return
end if
- if (Pcnst /= num_scalars) then
- write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, num_scalars
+ if (Pcnst /= block % state % time_levs(1) % state % num_scalars) then
+ write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, block % state % time_levs(1) % state % num_scalars
return
end if
+ uReconstZonal => block % diag % uReconstructZonal % array
+ uReconstMeridional => block % diag % uReconstructMeridional % array
+ scalars => block % state % time_levs(1) % state % scalars % array
+ theta => block % diag % theta % array
+ pressure_base => block % diag % pressure_base % array
+ pressure_p => block % diag % pressure_p % array
!
- ! Fill in CAM arrays from block % time_levs(1) % state arrays
+ ! Fill in CAM arrays from block % state % time_levs(1) % state arrays
!
do iCell=1,block % mesh % nCellsSolve
- Psd(iCell) = block % time_levs(1) % state % surface_pressure % array(iCell)
- Phis(iCell) = block % time_levs(1) % state % geopotential % array(1,iCell)
+ Psd(iCell) = block % diag % surface_pressure % array(iCell)
+ Phis(iCell) = block % mesh % zgrid % array(1,iCell) * 9.81
do k=1,block % mesh % nVertLevels
- T(iCell,k) = theta(k,iCell) * (0.5*(pressure(k,iCell)+pressure(k+1,iCell)) / p0) ** (rgas/cp)
- Omega(iCell,k) = 0.5*(ww(k,iCell) + ww(k+1,iCell))
+ theta(k,iCell) = block % state % time_levs(1) % state % theta_m % array(k,iCell) &
+ / (1.0 + 1.61 * scalars(index_qv,k,iCell))
+ T(iCell,k) = theta(k,iCell) * ((pressure_base(k,iCell)+pressure_p(k,iCell)) / p0) ** (rgas/cp)
+ Pmid(iCell,k) = pressure_base(k,iCell)+pressure_p(k,iCell)
+ if (k == 1) then
+ Pint(iCell,k) = Psd(iCell)
+ else
+ Pint(iCell,k) = block % mesh % fzm % array(k) * Pmid(iCell,k-1) + &
+ block % mesh % fzp % array(k) * Pmid(iCell,k)
+ end if
+ Zint(iCell,k) = block % mesh % zgrid % array(k,iCell)
+ Zmid(iCell,k) = 0.5*(block % mesh % zgrid % array(k,iCell) + block % mesh % zgrid % array(k+1,iCell))
+
+! Omega(iCell,k) = 0.5*(ww(k,iCell) + ww(k+1,iCell))
+ Omega(iCell,k) = 0.0 !MGD NEEDS FIXING
!
! NOTE: Eventually, we need to ensure that the moisture variables we return to CAM
! are what CAM expects, rather than simply what we have in our scalars array
!
- do iScalar=1,num_scalars
+ do iScalar=1,block % state % time_levs(1) % state % num_scalars
Tracer(iCell,k,iScalar) = scalars(iScalar,k,iCell) / (1.0 + scalars(index_qv,k,iCell))
end do
end do
+ Zint(iCell,Plev+1) = block % mesh % zgrid % array(Plev+1,iCell)
+ Pint(iCell,Plev+1) = Pmid(iCell,Plev) + (Pmid(iCell,Plev) - Pint(iCell,Plev))
end do
do iCell=1,block % mesh % nCellsSolve
@@ -687,8 +716,21 @@
Uy(iCell,k) = uReconstMeridional(k,iCell)
end do
end do
-#endif
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX uReconstructZonal=', minval(Ux), maxval(Ux)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX uReconstructMeridional=', minval(Uy), maxval(Uy)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Psd=', minval(Psd), maxval(Psd)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Phis=', minval(Phis), maxval(Phis)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Pint=', minval(Pint), maxval(Pint)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Pmid=', minval(Pmid), maxval(Pmid)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Zint=', minval(Zint), maxval(Zint)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Zmid=', minval(Zmid), maxval(Zmid)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX T=', minval(T), maxval(T)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Omega=', minval(Omega), maxval(Omega)
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Tracer1=', minval(Tracer(:,:,1)), maxval(Tracer(:,:,1))
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Tracer2=', minval(Tracer(:,:,2)), maxval(Tracer(:,:,2))
+write(0,*) '<- MPAS DEBUGGING: MIN/MAX Tracer3=', minval(Tracer(:,:,3)), maxval(Tracer(:,:,3))
+
end subroutine mpas_to_cam
@@ -717,34 +759,23 @@
real (kind=RKIND), dimension(Numcols,Plev,Pcnst), intent(in) :: Tracer
integer :: iCell, iEdge, iScalar, k, j
+ integer :: index_qv
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: edgesOnCell
- real (kind=RKIND), dimension(:,:), pointer :: h
+ real (kind=RKIND), dimension(:,:), pointer :: pressure_base, pressure_p
real (kind=RKIND), dimension(:,:), pointer :: east, north, edge_normal, pressure
- real (kind=RKIND), dimension(:,:), pointer :: Ux_tend_halo, Uy_tend_halo, theta_tend, u_tend
- real (kind=RKIND), dimension(:,:), pointer :: Ux_phys_tend, Uy_phys_tend
+ real (kind=RKIND), dimension(:,:), pointer :: theta_tend, u_tend
+ real (kind=RKIND), dimension(:,:), pointer :: zonal_tend, meridional_tend
real (kind=RKIND), dimension(:,:,:), pointer :: scalars, scalars_tend
type (block_type), pointer :: block
write(0,*) 'Called CAM_TO_MPAS'
-#if 0
block => domain % blocklist
- nEdgesOnCell => block % mesh % nEdgesOnCell % array
- edgesOnCell => block % mesh % edgesOnCell % array
- east => block % mesh % east % array
- north => block % mesh % north % array
- edge_normal => block % mesh % edgeNormalVectors % array
- pressure => block % time_levs(1) % state % pressure % array
- scalars => block % time_levs(1) % state % scalars % array
- scalars_tend => block % time_levs(1) % state % scalars_phys_tend % array
- theta_tend => block % time_levs(1) % state % theta_phys_tend % array
- u_tend => block % time_levs(1) % state % u_phys_tend % array
- Ux_phys_tend => block % time_levs(1) % state % Ux_phys_tend % array
- Uy_phys_tend => block % time_levs(1) % state % Uy_phys_tend % array
- h => block % time_levs(1) % state % h % array
+ index_qv = block % state % time_levs(1) % state % index_qv
+
!
! Perform basic sanity check on expected and available field dimensions
!
@@ -756,74 +787,73 @@
write(0,*) 'Error: mismatch between Plev and nVertLevels: ', Plev, block % mesh % nVertLevels
return
end if
- if (Pcnst /= num_scalars) then
- write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, num_scalars
+ if (Pcnst /= block % state % time_levs(1) % state % num_scalars) then
+ write(0,*) 'Error: mismatch between Pcnst and num_scalars: ', Pcnst, block % state % time_levs(1) % state % num_scalars
return
end if
+ nEdgesOnCell => block % mesh % nEdgesOnCell % array
+ edgesOnCell => block % mesh % edgesOnCell % array
+ east => block % mesh % east % array
+ north => block % mesh % north % array
+ edge_normal => block % mesh % edgeNormalVectors % array
+ scalars => block % state % time_levs(1) % state % scalars % array
+ pressure_base => block % diag % pressure_base % array
+ pressure_p => block % diag % pressure_p % array
+ scalars_tend => block % tend_physics % scalars % array
+ theta_tend => block % tend_physics % theta % array
+ u_tend => block % tend_physics % u % array
+ zonal_tend => block % tend_physics % ux % array
+ meridional_tend => block % tend_physics % uy % array
+
!
! Fill in MPAS tendency arrays from arguments
!
- allocate(Ux_tend_halo(block % mesh % nVertLevels, block % mesh % nCells))
- allocate(Uy_tend_halo(block % mesh % nVertLevels, block % mesh % nCells))
-
do iCell=1,block % mesh % nCellsSolve
do k=1,block % mesh % nVertLevels
- Ux_tend_halo(k,iCell) = Ux_tend(iCell,k)
- Uy_tend_halo(k,iCell) = Uy_tend(iCell,k)
+ zonal_tend(k,iCell) = Ux_tend(iCell,k)
+ meridional_tend(k,iCell) = Uy_tend(iCell,k)
end do
end do
- call dmpar_exch_halo_field2dReal(domain % dminfo, Ux_tend_halo, &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, Uy_tend_halo, &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call mpas_dmpar_exch_halo_field(block % tend_physics % ux)
+ call mpas_dmpar_exch_halo_field(block % tend_physics % uy)
u_tend(:,:) = 0.0
do iCell=1,block % mesh % nCells
do j=1,nEdgesOnCell(iCell)
iEdge = edgesOnCell(j,iCell)
do k=1,block % mesh % nVertLevels
- u_tend(k,iEdge) = u_tend(k,iEdge) + 0.5 * Ux_tend_halo(k,iCell) * (edge_normal(1,iEdge) * east(1,iCell) + &
- edge_normal(2,iEdge) * east(2,iCell) + &
- edge_normal(3,iEdge) * east(3,iCell)) &
- + 0.5 * Uy_tend_halo(k,iCell) * (edge_normal(1,iEdge) * north(1,iCell) + &
- edge_normal(2,iEdge) * north(2,iCell) + &
- edge_normal(3,iEdge) * north(3,iCell))
+ u_tend(k,iEdge) = u_tend(k,iEdge) + 0.5 * zonal_tend(k,iCell) * (edge_normal(1,iEdge) * east(1,iCell) + &
+ edge_normal(2,iEdge) * east(2,iCell) + &
+ edge_normal(3,iEdge) * east(3,iCell)) &
+ + 0.5 * meridional_tend(k,iCell) * (edge_normal(1,iEdge) * north(1,iCell) + &
+ edge_normal(2,iEdge) * north(2,iCell) + &
+ edge_normal(3,iEdge) * north(3,iCell))
end do
end do
- ! Store original zonal and meridional tendencies from CAM to be written to output files for testing
- do k=1,block % mesh % nVertLevels
- Ux_phys_tend(k,iCell) = Ux_tend_halo(k,iCell)
- Uy_phys_tend(k,iCell) = Uy_tend_halo(k,iCell)
- end do
end do
- deallocate(Ux_tend_halo)
- deallocate(Uy_tend_halo)
-
do iCell=1,block % mesh % nCellsSolve
do k=1,block % mesh % nVertLevels
- theta_tend(k,iCell) = T_tend(iCell,k) * (p0 / (0.5*(pressure(k,iCell)+pressure(k+1,iCell)))) ** (rgas/cp)
+ theta_tend(k,iCell) = T_tend(iCell,k) * (p0 / (pressure_base(k,iCell)+pressure_p(k,iCell))) ** (rgas/cp)
! Couple theta tendency with h valid at time when tendency was computed
- theta_tend(k,iCell) = theta_tend(k,iCell) * h(k,iCell)
+! theta_tend(k,iCell) = theta_tend(k,iCell) * h(k,iCell)
!
! NOTE: Once we begin to use more than just qv in MPAS, we need to make sure to
! properly assign variables from Tracer array provided from CAM physics
!
- do iScalar=1,num_scalars
+ do iScalar=1,block % state % time_levs(1) % state % num_scalars
scalars_tend(iScalar,k,iCell) = (Tracer(iCell,k,iScalar) * (1.0 + Tracer(iCell,k,index_qv)) &
- scalars(iScalar,k,iCell)&
) / dt_physics
! Couple scalar tendencies with h valid at time when tendencies were computed
- scalars_tend(iScalar,k,iCell) = scalars_tend(iScalar,k,iCell) * h(k,iCell)
+! scalars_tend(iScalar,k,iCell) = scalars_tend(iScalar,k,iCell) * h(k,iCell)
end do
end do
end do
@@ -832,16 +862,15 @@
!
! Do a ghost-cell update on tendency arrays
!
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % theta_phys_tend % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % u_phys_tend % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(1) % state % scalars_phys_tend % array(:,:,:), &
- num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-#endif
+ call mpas_dmpar_exch_halo_field(block % tend_physics % theta)
+ call mpas_dmpar_exch_halo_field(block % tend_physics % u)
+ call mpas_dmpar_exch_halo_field(block % tend_physics % scalars)
+
+write(0,*) '-> MPAS DEBUGGING: MIN/MAX u_tend=', minval(u_tend(:,1:block % mesh % nEdges)),maxval(u_tend(:,1:block % mesh % nEdges))
+write(0,*) '-> MPAS DEBUGGING: MIN/MAX theta_tend=', minval(theta_tend(:,1:block % mesh % nCells)),maxval(theta_tend(:,1:block % mesh % nCells))
+write(0,*) '-> MPAS DEBUGGING: MIN/MAX qv_tend=', minval(scalars_tend(1,:,1:block % mesh % nCells)),maxval(scalars_tend(1,:,1:block % mesh % nCells))
+write(0,*) '-> MPAS DEBUGGING: MIN/MAX qc_tend=', minval(scalars_tend(2,:,1:block % mesh % nCells)),maxval(scalars_tend(2,:,1:block % mesh % nCells))
+write(0,*) '-> MPAS DEBUGGING: MIN/MAX qr_tend=', minval(scalars_tend(3,:,1:block % mesh % nCells)),maxval(scalars_tend(3,:,1:block % mesh % nCells))
end subroutine cam_to_mpas
@@ -855,25 +884,80 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
-
+
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block_ptr
+
+ type (MPAS_Time_Type) :: currTime
+ character(len=StrKIND) :: timeStamp
+ integer :: itimestep
+ integer :: ierr
integer :: idynstep
write(0,*) 'Called MPAS_DYN_RUN'
-#if 0
+
+ ! Eventually, dt should be domain specific
+ dt = config_dt
+
+ ! During integration, time level 1 stores the model state at the beginning of the
+ ! time step, and time level 2 stores the state advanced dt_dynamics in time by timestep(...)
+write(0,*) 'MGD running dynamics for n_subcycle_steps=',n_subcycle_steps
do idynstep=1,n_subcycle_steps
+
+ currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
+ call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
+ write(0,*) 'Begin timestep ', trim(timeStamp)
+
+ ! Input external updates (i.e. surface)
+ if (mpas_is_alarm_ringing(clock, sfcAlarmID, ierr=ierr)) then
+ call mpas_reset_clock_alarm(clock, sfcAlarmID, ierr=ierr)
+
+ call mpas_read_and_distribute_fields(sfc_update_obj)
+ sfc_update_obj % time = sfc_update_obj % time + 1
+ end if
+
+ call mpas_timer_start("time integration")
+ call atm_do_timestep(domain, dt_dynamics, itimestep)
+ call mpas_timer_stop("time integration")
+
+ ! Move time level 2 fields back into time level 1 for next time step
+ call mpas_shift_time_levels_state(domain % blocklist % state)
+
+
+ ! Advance clock before writing output
itimestep = itimestep + 1
- call mpas_timestep(domain, itimestep, dt_dynamics)
+ call mpas_advance_clock(clock)
+ currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
+ call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
+
+ !TODO: MPAS_getClockRingingAlarms is probably faster than multiple MPAS_isAlarmRinging...
+ if (mpas_is_alarm_ringing(clock, outputAlarmID, ierr=ierr)) then
+ call mpas_reset_clock_alarm(clock, outputAlarmID, ierr=ierr)
+ ! output_frame will always be > 1 here unless it was reset after the maximum number of frames per outfile was reached
+ if(output_frame == 1) then
+ call mpas_output_state_finalize(output_obj, domain % dminfo)
+ call mpas_output_state_init(output_obj, domain, "OUTPUT", trim(timeStamp))
+ end if
+ call atm_write_output_frame(output_obj, output_frame, domain)
+ end if
+
+ if (mpas_is_alarm_ringing(clock, restartAlarmID, ierr=ierr)) then
+ call mpas_reset_clock_alarm(clock, restartAlarmID, ierr=ierr)
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call atm_compute_restart_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % diag, block_ptr % mesh)
+ block_ptr => block_ptr % next
+ end do
- ! Move time level 2 fields back into time level 1 for next time step
- call shift_time_levels(domain)
+ ! Write one restart time per file
+ call mpas_output_state_init(restart_obj, domain, "RESTART", trim(timeStamp))
+ call mpas_output_state_for_domain(restart_obj, domain, 1)
+ call mpas_output_state_finalize(restart_obj, domain % dminfo)
+ end if
- if (mod(itimestep, config_output_interval) == 0) then
- output_frame = output_frame + 1
- call output_state_for_domain(output_obj, domain, output_frame)
- end if
end do
-#endif
end subroutine mpas_dyn_run
@@ -891,12 +975,12 @@
write(0,*) 'Called MPAS_FINAL'
-#if 0
- call output_state_finalize(output_obj, domain % dminfo)
- call deallocate_domain(domain)
+ call mpas_output_state_finalize(output_obj, domain % dminfo)
- call dmpar_finalize(dminfo)
-#endif
+ call mpas_core_finalize(domain)
+
+ call mpas_framework_finalize(dminfo, domain)
+
end subroutine mpas_final
</font>
</pre>