<p><b>duda</b> 2012-05-30 12:32:23 -0600 (Wed, 30 May 2012)</p><p>BRANCH COMMIT<br>
<br>
Commit placeholder interface module used by CAM to interact with MPAS-A.<br>
<br>
<br>
A src/driver/mpas_cam_interface.F<br>
M src/driver/Makefile<br>
</p><hr noshade><pre><font color="gray">Modified: branches/cam_mpas_nh/src/driver/Makefile
===================================================================
--- branches/cam_mpas_nh/src/driver/Makefile 2012-05-30 17:21:57 UTC (rev 1947)
+++ branches/cam_mpas_nh/src/driver/Makefile 2012-05-30 18:32:23 UTC (rev 1948)
@@ -1,12 +1,15 @@
.SUFFIXES: .F .o
OBJS = mpas_subdriver.o \
+ mpas_cam_interface.o \
mpas.o
all: $(OBJS)
mpas_subdriver.o:
+mpas_cam_interface.o: mpas_subdriver.o
+
mpas.o: mpas_subdriver.o
clean:
Added: branches/cam_mpas_nh/src/driver/mpas_cam_interface.F
===================================================================
--- branches/cam_mpas_nh/src/driver/mpas_cam_interface.F (rev 0)
+++ branches/cam_mpas_nh/src/driver/mpas_cam_interface.F 2012-05-30 18:32:23 UTC (rev 1948)
@@ -0,0 +1,976 @@
+module mpas_cam_interface
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_io_input
+ use mpas_io_output
+ use mpas_dmpar
+ use mpas_constants
+
+
+ type (dm_info), pointer :: dminfo
+ type (domain_type), pointer :: domain
+ real (kind=RKIND) :: dt_dynamics, dt_physics, p0
+ integer :: itimestep, n_subcycle_steps
+ type (io_output_object), save :: output_obj
+ logical :: restart_run
+
+
+ contains
+
+
+ subroutine mpas_init1(mpi_comm, phys_dt)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE MPAS_INIT1
+ !
+ ! Initializes the MPAS software infrastructure, reads run-time configuration
+ ! information, reads grid information from an MPAS grid file, and allocates
+ ! storage for fields to be provided by CAM through either
+ ! cam_inidat_to_mpas() or cam_restart_to_mpas().
+ !
+ ! Input: mpi_comm - an MPI communicator supplied by CAM to be used by MPAS
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer, intent(in) :: mpi_comm
+ real (kind=RKIND), intent(in) :: phys_dt
+
+ integer :: iCell, iEdge, j
+ real (kind=RKIND), dimension(:), pointer :: latCell, lonCell
+ real (kind=RKIND), dimension(:,:), pointer :: east, north
+
+ write(0,*) 'Called MPAS_INIT1'
+
+#if 0
+ allocate(dminfo)
+ call dmpar_init(dminfo, mpi_comm)
+
+ call read_namelist(dminfo)
+
+ itimestep = 0
+ restart_run = .false. ! If a restart run, this will be set to true in cam_restart_to_mpas
+
+ !
+ ! Set physics timestep and verify that it is evenly divided by dynamics timestep
+ !
+ dt_physics = phys_dt
+ dt_dynamics = dt_physics / ceiling(dt_physics / config_dt)
+ 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)
+
+ latCell => domain % blocklist % mesh % latCell % array
+ lonCell => domain % blocklist % mesh % lonCell % array
+ east => domain % blocklist % mesh % east % array
+ north => domain % blocklist % mesh % north % array
+
+
+ !
+ ! Compute unit vectors in east and north directions for each cell
+ !
+ do iCell = 1,domain % blocklist % mesh % nCellsSolve
+
+ east(1,iCell) = -sin(lonCell(iCell))
+ east(2,iCell) = cos(lonCell(iCell))
+ east(3,iCell) = 0.0
+ call r3_normalize(east(1,iCell), east(2,iCell), east(3,iCell))
+
+ north(1,iCell) = -cos(lonCell(iCell))*sin(latCell(iCell))
+ north(2,iCell) = -sin(lonCell(iCell))*sin(latCell(iCell))
+ north(3,iCell) = cos(latCell(iCell))
+ call r3_normalize(north(1,iCell), north(2,iCell), north(3,iCell))
+
+ end do
+#endif
+
+ end subroutine mpas_init1
+
+
+ subroutine cam_inidat_to_mpas(Numcols, MaxEdges, Plev, Pcnst, Pref, Ptop, &
+ Ae, Be, Ac, Bc, Psd, Phis, Theta, U, Tracer)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE CAM_INIDAT_TO_MPAS
+ !
+ ! Input: Numcols - number of columns/cells
+ ! MaxEdges - maximum number of edges per cell
+ ! Plev - number of vertical levels
+ ! Pcnst - number of tracers
+ ! Pref - reference surface pressure (Pa)
+ ! Ptop - pressure at top of model (Pa)
+ ! Ae, Be - A and B arrays at vertical edges
+ ! Ac, Bc - A and B arrays at vertical centers
+ ! Psd - dry surface pressure (Pa)
+ ! Phis - surface geopotential (m^2/s^2)
+ ! Theta - potential temperature (K)
+ ! U - normal velocity at edges (m/s)
+ ! Tracer - tracer mixing ratios (kg/kg)
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer, intent(in) :: Numcols
+ integer, intent(in) :: MaxEdges
+ integer, intent(in) :: Plev
+ integer, intent(in) :: Pcnst
+ real (kind=RKIND), intent(in) :: Pref
+ real (kind=RKIND), intent(in) :: Ptop
+ real (kind=RKIND), dimension(Plev+1), intent(in) :: Ae, Be
+ real (kind=RKIND), dimension(Plev), intent(in) :: Ac, Bc
+ real (kind=RKIND), dimension(Numcols), intent(in) :: Psd
+ real (kind=RKIND), dimension(Numcols), intent(in) :: Phis
+ real (kind=RKIND), dimension(Numcols,Plev), intent(in) :: Theta
+ real (kind=RKIND), dimension(Numcols,MaxEdges,Plev), intent(in) :: U
+ real (kind=RKIND), dimension(Numcols,Plev,Pcnst), intent(in) :: Tracer
+
+ integer :: iCell, iEdge, iScalar, k, i
+ real (kind=RKIND) :: qtot
+ real (kind=RKIND), dimension(Plev+1) :: znu, znw, znwc, znwv
+ real (kind=RKIND), dimension(Plev+1) :: znuc, znuv, bn, divh, dpn
+ type (block_type), pointer :: block
+ real (kind=RKIND), dimension(:), pointer :: surface_pressure, rdnu, rdnw, fnm, fnp, dbn, dnu, dnw
+ real (kind=RKIND), dimension(:,:), pointer :: t_pot, u_normal, pressure, geopotential, alpha, h
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars
+
+ write(0,*) 'Called CAM_INIDAT_TO_MPAS'
+
+#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
+ rdnu => block % mesh % rdnu % array
+ rdnw => block % mesh % rdnw % array
+ fnm => block % mesh % fnm % array
+ fnp => block % mesh % fnp % array
+ dbn => block % mesh % dbn % array
+ dnu => block % mesh % dnu % array
+ dnw => block % mesh % dnw % array
+
+
+ !
+ ! Perform basic sanity check on expected and available field dimensions
+ !
+ if (Numcols /= block % mesh % nCellsSolve) then
+ write(0,*) 'Error: mismatch between Numcols and nCellsSolve: ', Numcols, block % mesh % nCellsSolve
+ return
+ end if
+ if (Plev /= block % mesh % nVertLevels) then
+ 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
+ return
+ end if
+
+ write(0,*) 'Numcols = ', Numcols
+ write(0,*) 'MaxEdges = ', MaxEdges
+ write(0,*) 'Plev = ', Plev
+ write(0,*) 'Pcnst = ', Pcnst
+ write(0,*) 'Pref = ', Pref
+ write(0,*) 'Ptop = ', Ptop
+ write(0,*) 'min/max for Ae = ', minval(Ae), maxval(Ae)
+ write(0,*) 'min/max for Be = ', minval(Be), maxval(Be)
+ write(0,*) 'min/max for Ac = ', minval(Ac), maxval(Ac)
+ write(0,*) 'min/max for Bc = ', minval(Bc), maxval(Bc)
+ write(0,*) 'min/max for Psd = ', minval(Psd), maxval(Psd)
+ write(0,*) 'min/max for Phis = ', minval(Phis), maxval(Phis)
+ write(0,*) 'min/max for Theta = ', minval(Theta), maxval(Theta)
+ write(0,*) 'min/max for U = ', minval(U), maxval(U)
+ write(0,*) 'min/max for Tracer = ', minval(Tracer), maxval(Tracer)
+
+
+ !
+ ! Scale geometry information from unit sphere to sphere with radius a (Earth radius defined in constants module)
+ !
+ block % mesh % xCell % array = block % mesh % xCell % array * a
+ block % mesh % yCell % array = block % mesh % yCell % array * a
+ block % mesh % zCell % array = block % mesh % zCell % array * a
+ block % mesh % xVertex % array = block % mesh % xVertex % array * a
+ block % mesh % yVertex % array = block % mesh % yVertex % array * a
+ block % mesh % zVertex % array = block % mesh % zVertex % array * a
+ block % mesh % xEdge % array = block % mesh % xEdge % array * a
+ block % mesh % yEdge % array = block % mesh % yEdge % array * a
+ block % mesh % zEdge % array = block % mesh % zEdge % array * a
+ block % mesh % dvEdge % array = block % mesh % dvEdge % array * a
+ block % mesh % dcEdge % array = block % mesh % dcEdge % array * a
+ block % mesh % areaCell % array = block % mesh % areaCell % array * a**2.0
+ block % mesh % areaTriangle % array = block % mesh % areaTriangle % array * a**2.0
+ block % mesh % kiteAreasOnVertex % array = block % mesh % kiteAreasOnVertex % array * a**2.0
+
+
+ !
+ ! Derive vertical coordinate information
+ !
+ p0 = Pref
+ bn (1) = 1.
+ znw(1) = 1.
+ znwc(1) = 1.
+ znwv(1) = (znwc(1)-.252)*pii/2.
+
+ do k=1,Plev
+ znuc(Plev+1-k) = Ac(Plev+1-k) + Bc(Plev+1-k)
+ znwc(Plev+1-k+1) = Ae(Plev+2-k) + Be(Plev+2-k)
+ znu (Plev+1-k ) = (znuc(Plev+1-k )*p0-Ptop)/(p0-Ptop)
+ znw (Plev+1-k+1) = (znwc(Plev+1-k+1)*p0-Ptop)/(p0-Ptop)
+ bn(k+1) = Be(k+1)
+ end do
+
+ do k=1,Plev
+ znuv(k ) = (znuc(k )-.252)*pii/2.
+ znwv(k+1) = (znwc(k+1)-.252)*pii/2.
+ dnw (k) = znw(k+1)-znw(k)
+ rdnw(k) = 1./dnw(k)
+ dbn (k) = rdnw(k)*(bn(k+1)-bn(k))
+ dpn (k) = 0.
+ divh(k) = 0.
+ end do
+
+ dpn(Plev+1)=0.
+ fnm(1) = 0.
+ fnp(1) = 0.
+ do k=2,Plev
+ dnu (k) = .5*(dnw(k)+dnw(k-1))
+ rdnu(k) = 1./dnu(k)
+ fnp (k) = .5* dnw(k )/dnu(k)
+ fnm (k) = .5* dnw(k-1)/dnu(k)
+ end do
+
+ !
+ ! Initialize u_normal
+ !
+ u_normal(:,:) = 0.0
+ do iCell=1,Numcols
+ do i=1,block % mesh % nEdgesOnCell % array(iCell)
+ iEdge = block % mesh % edgesOnCell % array(i,iCell)
+ do k=1,Plev
+ u_normal(k,iEdge) = U(iCell,i,k)
+ end do
+ end do
+ end do
+
+ !
+ ! Initialize surface pressure
+ !
+ surface_pressure(:) = 0.0
+ do iCell=1,Numcols
+ surface_pressure(iCell) = Psd(iCell)
+ end do
+
+ !
+ ! Initialize pressure, theta, alpha, and geopotential
+ !
+ h(:,:) = 0.0
+ pressure(:,:) = 0.0
+ t_pot(:,:) = 0.0
+ alpha(:,:) = 0.0
+ geopotential(:,:) = 0.0
+ do iCell=1,Numcols
+ do k=1,Plev
+ h(k,iCell) = (1.-dbn(k))*(p0-Ptop)+dbn(k)*(surface_pressure(iCell)-Ptop)
+ end do
+
+ pressure(Plev+1,iCell) = Ptop
+ do k=Plev,1,-1
+ qtot = sum(Tracer(iCell,k,1:Pcnst))
+ pressure(k,iCell) = pressure(k+1,iCell)-dnw(k)*h(k,iCell)*(1.0 + qtot)
+ end do
+
+ do k=1,Plev
+ t_pot(k,iCell) = Theta(iCell,k)
+ alpha(k,iCell) = (rgas/p0) * Theta(iCell,k) * (1.0 + 1.61*Tracer(iCell,k,1)) & ! NB: Assume that first Tracer is qv
+ * (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**cvpm
+ end do
+
+ geopotential(1,iCell) = Phis(iCell)
+ do k=1,Plev
+ geopotential(k+1,iCell) = geopotential(k,iCell)-dnw(k)*h(k,iCell)*alpha(k,iCell)
+ end do
+ end do
+
+ !
+ ! Initialize scalars
+ !
+ scalars(:,:,:) = 0.0
+ do iCell=1,Numcols
+ do k=1,Plev
+ do iScalar=1,Pcnst
+ scalars(iScalar,k,iCell) = Tracer(iCell,k,iScalar)
+ end do
+ end do
+ end do
+
+
+ !
+ ! 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)
+
+ do i=2,nTimeLevs
+ call copy_state(block % time_levs(1) % state, block % time_levs(i) % state)
+ end do
+ block % time_levs(1) % state % xtime % scalar = 0.0
+#endif
+
+ end subroutine cam_inidat_to_mpas
+
+
+ subroutine mpas_restart_to_cam(Numcols, MaxEdges, Plev, Pcnst, Psd, Phis, Theta, U, W, Tracer)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE MPAS_RESTART_TO_CAM
+ !
+ ! Input: Numcols - number of columns/cells
+ ! MaxEdges - maximum number of edges per cell
+ ! Plev - number of vertical levels
+ ! Pcnst - number of tracers
+ ! Psd - dry surface pressure (Pa)
+ ! Phis - surface geopotential (m^2/s^2)
+ ! Theta - potential temperature (K)
+ ! U - normal velocity at edges (m/s)
+ ! W - vertical velocity (Pa/s)
+ ! Tracer - tracer mixing ratios (kg/kg)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer, intent(in) :: Numcols
+ integer, intent(in) :: MaxEdges
+ integer, intent(in) :: Plev
+ 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), intent(out) :: Theta
+ real (kind=RKIND), dimension(Numcols,MaxEdges,Plev), intent(out) :: U
+ real (kind=RKIND), dimension(Numcols,Plev+1), intent(out) :: W
+ real (kind=RKIND), dimension(Numcols,Plev,Pcnst), intent(out) :: Tracer
+
+ integer :: iCell, iEdge, k, iScalar, i
+ type (block_type), pointer :: block
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnCell
+
+ write(0,*) 'Called MPAS_RESTART_TO_CAM'
+
+#if 0
+ block => domain % blocklist
+
+ nEdgesOnCell => block % mesh % nEdgesOnCell % array
+ edgesOnCell => block % mesh % edgesOnCell % array
+
+ !
+ ! Perform basic sanity check on expected and available field dimensions
+ !
+ if (Numcols /= block % mesh % nCellsSolve) then
+ write(0,*) 'Error: mismatch between Numcols and nCellsSolve: ', Numcols, block % mesh % nCellsSolve
+ return
+ end if
+ if (Plev /= block % mesh % nVertLevels) then
+ 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
+ return
+ end if
+
+
+ !
+ ! Fill in output arrays from block % 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)
+ 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)
+ !
+ ! 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)
+ end do
+ do i=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i,iCell)
+ U(iCell,i,k) = block % time_levs(1) % state % u % array(k,iEdge)
+ end do
+ end do
+ W(iCell,k) = block % time_levs(1) % state % ww % array(k,iCell)
+ end do
+#endif
+
+ end subroutine mpas_restart_to_cam
+
+
+ subroutine cam_restart_to_mpas(Numcols, MaxEdges, Plev, Pcnst, Pref, Ptop, &
+ Ae, Be, Ac, Bc, Psd, Phis, Theta, U, W, Tracer)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE CAM_RESTART_TO_MPAS
+ !
+ ! Input: Numcols - number of columns/cells
+ ! MaxEdges - maximum number of edges per cell
+ ! Plev - number of vertical levels
+ ! Pcnst - number of tracers
+ ! Pref - reference surface pressure (Pa)
+ ! Ptop - pressure at top of model (Pa)
+ ! Ae, Be - A and B arrays at vertical edges
+ ! Ac, Bc - A and B arrays at vertical centers
+ ! Psd - dry surface pressure (Pa)
+ ! Phis - surface geopotential (m^2/s^2)
+ ! Theta - potential temperature (K)
+ ! U - normal velocity at edges (m/s)
+ ! W - vertical velocity (Pa/s)
+ ! Tracer - tracer mixing ratios (kg/kg)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer, intent(in) :: Numcols
+ integer, intent(in) :: MaxEdges
+ integer, intent(in) :: Plev
+ integer, intent(in) :: Pcnst
+ real (kind=RKIND), intent(in) :: Pref
+ real (kind=RKIND), intent(in) :: Ptop
+ real (kind=RKIND), dimension(Plev+1), intent(in) :: Ae, Be
+ real (kind=RKIND), dimension(Plev), intent(in) :: Ac, Bc
+ real (kind=RKIND), dimension(Numcols), intent(in) :: Psd
+ real (kind=RKIND), dimension(Numcols), intent(in) :: Phis
+ real (kind=RKIND), dimension(Numcols,Plev), intent(in) :: Theta
+ real (kind=RKIND), dimension(Numcols,MaxEdges,Plev), intent(in) :: U
+ real (kind=RKIND), dimension(Numcols,Plev+1), intent(in) :: W
+ real (kind=RKIND), dimension(Numcols,Plev,Pcnst), intent(in) :: Tracer
+
+ integer :: iCell, k, i
+ type (block_type), pointer :: block
+
+ restart_run = .true.
+
+ write(0,*) 'Called CAM_RESTART_TO_MPAS'
+
+#if 0
+ block => domain % blocklist
+
+ !
+ ! Perform basic sanity check on expected and available field dimensions
+ !
+ if (Numcols /= block % mesh % nCellsSolve) then
+ write(0,*) 'Error: mismatch between Numcols and nCellsSolve: ', Numcols, block % mesh % nCellsSolve
+ return
+ end if
+ if (Plev /= block % mesh % nVertLevels) then
+ 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
+ return
+ end if
+
+
+ !
+ ! Fill in MPAS fields in block % 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, &
+ Ae, Be, Ac, Bc, Psd, Phis, Theta, U, Tracer)
+
+ do iCell=1,Numcols
+ do k=1,Plev+1
+ block % time_levs(1) % state % ww % array(k,iCell) = W(iCell,k)
+ end do
+ end do
+
+
+ !
+ ! Do a ghost-cell update
+ !
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % 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(:,:)
+ end do
+#endif
+
+ end subroutine cam_restart_to_mpas
+
+
+ subroutine mpas_init2()
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE MPAS_INIT2
+ !
+ ! Calls the core-specific (most likely the hydrostatic atmosphere core)
+ ! initialization routine mpas_init() after initial fields have been provided
+ ! through a call to either cam_inidat_to_mpas() or cam_restart_to_mpas().
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer :: i, iCell, k, iScalar
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars
+ type (block_type), pointer :: block
+
+ write(0,*) 'Called MPAS_INIT2'
+
+#if 0
+ !
+ ! If this is not a restart run, we first need to transform the initial scalar fields
+ ! from wet mixing ratios to dry mixing ratios
+ !
+ if (.not. restart_run) then
+ block => domain % blocklist
+ do while (associated(block))
+ scalars => block % time_levs(1) % state % scalars % array
+
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ do iScalar=1,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(:,:,:)
+ end do
+
+ block => block % next
+ end do
+ end if
+
+ ! 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)
+ 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
+
+ end subroutine mpas_init2
+
+
+ subroutine mpas_to_cam(Numcols, Plev, Pcnst, Psd, Phis, T, Ux, Uy, Omega, Tracer)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE MPAS_TO_CAM
+ !
+ ! Input: Numcols - number of columns/cells
+ ! Plev - number of vertical levels
+ ! Pcnst - number of tracers
+ !
+ ! Output: Psd - dry surface pressure (Pa)
+ ! Phis - surface geopotential (m^2/s^2)
+ ! T - temperature (K)
+ ! Ux - longitudinal velocity at cell centers (m/s)
+ ! Uy - latitudinal velocity at cell centers (m/s)
+ ! Omega - omega (Pa/s)
+ ! Tracer - tracer mixing ratios (kg/kg)
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer, intent(in) :: Numcols
+ integer, intent(in) :: Plev
+ 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), intent(out) :: T
+ real (kind=RKIND), dimension(Numcols,Plev), intent(out) :: Ux
+ real (kind=RKIND), dimension(Numcols,Plev), intent(out) :: Uy
+ real (kind=RKIND), dimension(Numcols,Plev), intent(out) :: Omega
+ real (kind=RKIND), dimension(Numcols,Plev,Pcnst), intent(out) :: Tracer
+
+ integer :: iCell, k, iScalar
+ real (kind=RKIND), dimension(:), pointer :: latCell, lonCell
+ real (kind=RKIND), dimension(:,:), pointer :: theta, pressure, 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
+
+
+ !
+ ! Perform basic sanity check on expected and available field dimensions
+ !
+ if (Numcols /= block % mesh % nCellsSolve) then
+ write(0,*) 'Error: mismatch between Numcols and nCellsSolve: ', Numcols, block % mesh % nCellsSolve
+ return
+ end if
+ if (Plev /= block % mesh % nVertLevels) then
+ 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
+ return
+ end if
+
+
+ !
+ ! Fill in CAM arrays from block % 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)
+ 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))
+ !
+ ! 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) = scalars(iScalar,k,iCell) / (1.0 + scalars(index_qv,k,iCell))
+ end do
+ end do
+ end do
+
+ do iCell=1,block % mesh % nCellsSolve
+ do k=1,block % mesh % nVertLevels
+ Ux(iCell,k) = uReconstZonal(k,iCell)
+ Uy(iCell,k) = uReconstMeridional(k,iCell)
+ end do
+ end do
+#endif
+
+ end subroutine mpas_to_cam
+
+
+ subroutine cam_to_mpas(Numcols, Plev, Pcnst, T_tend, Ux_tend, Uy_tend, Tracer)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE CAM_TO_MPAS
+ !
+ ! Input: Numcols - number of columns/cells
+ ! Plev - number of vertical levels
+ ! Pcnst - number of tracers
+ ! T_tend - temperature tendency (K/s)
+ ! Ux_tend - longitudinal velocity tendency at cell centers (m/s^2)
+ ! Uy_tend - latitudinal velocity tendency at cell centers (m/s^2)
+ ! Tracer - tracer mixing ratios (kg/kg)
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer, intent(in) :: Numcols
+ integer, intent(in) :: Plev
+ integer, intent(in) :: Pcnst
+ real (kind=RKIND), dimension(Numcols,Plev), intent(in) :: T_tend
+ real (kind=RKIND), dimension(Numcols,Plev), intent(in) :: Ux_tend
+ real (kind=RKIND), dimension(Numcols,Plev), intent(in) :: Uy_tend
+ real (kind=RKIND), dimension(Numcols,Plev,Pcnst), intent(in) :: Tracer
+
+ integer :: iCell, iEdge, iScalar, k, j
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnCell
+ real (kind=RKIND), dimension(:,:), pointer :: h
+ 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 :: 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
+
+ !
+ ! Perform basic sanity check on expected and available field dimensions
+ !
+ if (Numcols /= block % mesh % nCellsSolve) then
+ write(0,*) 'Error: mismatch between Numcols and nCellsSolve: ', Numcols, block % mesh % nCellsSolve
+ return
+ end if
+ if (Plev /= block % mesh % nVertLevels) then
+ 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
+ return
+ end if
+
+
+ !
+ ! 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)
+ 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)
+
+ 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))
+ 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)
+
+ ! Couple theta tendency with h valid at time when tendency was computed
+ 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
+ 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)
+ end do
+ end do
+ end do
+
+
+ !
+ ! 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
+
+ end subroutine cam_to_mpas
+
+
+ subroutine mpas_dyn_run()
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE MPAS_DYN_RUN
+ !
+ ! Advance MPAS dynamics by one timestep.
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer :: idynstep
+
+ write(0,*) 'Called MPAS_DYN_RUN'
+
+#if 0
+ do idynstep=1,n_subcycle_steps
+ itimestep = itimestep + 1
+ call mpas_timestep(domain, itimestep, dt_dynamics)
+
+ ! Move time level 2 fields back into time level 1 for next time step
+ call shift_time_levels(domain)
+
+ 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
+
+
+ subroutine mpas_final()
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE MPAS_FINAL
+ !
+ ! Deallocates storage allocated by MPAS and shuts down the MPAS software
+ ! infrastructure.
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ write(0,*) 'Called MPAS_FINAL'
+
+#if 0
+ call output_state_finalize(output_obj, domain % dminfo)
+ call deallocate_domain(domain)
+
+ call dmpar_finalize(dminfo)
+#endif
+
+ end subroutine mpas_final
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE ROTATE_ABOUT_VECTOR
+ !
+ ! Rotates the point (x,y,z) through an angle theta about the vector
+ ! originating at (a, b, c) and having direction (u, v, w).
+ !
+ ! Reference: http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/ArbitraryAxisRotation.html
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine rotate_about_vector(x, y, z, theta, a, b, c, u, v, w, xp, yp, zp)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: x, y, z, theta, a, b, c, u, v, w
+ real (kind=RKIND), intent(out) :: xp, yp, zp
+
+ real (kind=RKIND) :: vw2, uw2, uv2
+ real (kind=RKIND) :: m
+
+ vw2 = v**2.0 + w**2.0
+ uw2 = u**2.0 + w**2.0
+ uv2 = u**2.0 + v**2.0
+ m = sqrt(u**2.0 + v**2.0 + w**2.0)
+
+ xp = (a*vw2 + u*(-b*v-c*w+u*x+v*y+w*z) + ((x-a)*vw2+u*(b*v+c*w-v*y-w*z))*cos(theta) + m*(-c*v+b*w-w*y+v*z)*sin(theta))/m**2.0
+ yp = (b*uw2 + v*(-a*u-c*w+u*x+v*y+w*z) + ((y-b)*uw2+v*(a*u+c*w-u*x-w*z))*cos(theta) + m*( c*u-a*w+w*x-u*z)*sin(theta))/m**2.0
+ zp = (c*uv2 + w*(-a*u-b*v+u*x+v*y+w*z) + ((z-c)*uv2+w*(a*u+b*v-u*x-v*y))*cos(theta) + m*(-b*u+a*v-v*x+u*y)*sin(theta))/m**2.0
+
+ end subroutine rotate_about_vector
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE R3_CROSS
+ !
+ ! Computes the cross product of (ax, ay, az) and (bx, by, bz).
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine r3_cross(ax, ay, az, bx, by, bz, cx, cy, cz)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az
+ real (kind=RKIND), intent(in) :: bx, by, bz
+ real (kind=RKIND), intent(out) :: cx, cy, cz
+
+ cx = ay * bz - az * by
+ cy = az * bx - ax * bz
+ cz = ax * by - ay * bx
+
+ end subroutine r3_cross
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE R3_NORMALIZE
+ !
+ ! Normalizes the vector (ax, ay, az)
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine r3_normalize(ax, ay, az)
+
+ implicit none
+
+ real (kind=RKIND), intent(inout) :: ax, ay, az
+
+ real (kind=RKIND) :: mi
+
+ mi = 1.0 / sqrt(ax**2 + ay**2 + az**2)
+ ax = ax * mi
+ ay = ay * mi
+ az = az * mi
+
+ end subroutine r3_normalize
+
+end module mpas_cam_interface
</font>
</pre>