<p><b>duda</b> 2010-03-31 16:55:43 -0600 (Wed, 31 Mar 2010)</p><p>BRANCH COMMIT<br>
<br>
module_mpas_cam_interface.F:<br>
- Add implementation of CAM_INIDAT_TO_MPAS(...)<br>
<br>
- Add routine mpas_write_output() to allow us to <br>
write output.nc file from CAM test driver<br>
<br>
- Fix MPAS_DYN_RUN: after taking a timestep, we<br>
need to shift fields in time_levs(2) to<br>
time_levs(1)<br>
<br>
test_driver.F:<br>
- Rename variable 'Omega' to 'dPdt' to avoid name<br>
conflict with variable omega that gets use-associated<br>
from mpas_cam_interface via constants module<br>
<br>
- Run MPAS for 24 timesteps after calling CAM_INIDAT_TO_MPAS(...)<br>
and write an output file<br>
<br>
<br>
M src/driver_cam_interface/module_mpas_cam_interface.F<br>
M src/driver_cam_interface/test_driver.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/mpas_cam_coupling/src/driver_cam_interface/module_mpas_cam_interface.F
===================================================================
--- branches/mpas_cam_coupling/src/driver_cam_interface/module_mpas_cam_interface.F 2010-03-31 20:43:23 UTC (rev 175)
+++ branches/mpas_cam_coupling/src/driver_cam_interface/module_mpas_cam_interface.F 2010-03-31 22:55:43 UTC (rev 176)
@@ -5,6 +5,7 @@
use io_input
use dmpar
use subdriver
+ use constants
type (dm_info), pointer :: dminfo
@@ -97,36 +98,189 @@
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) :: p0
+ 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 :: theta, u_normal, pressure, geopotential, alpha, h
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars
write(0,*) 'Called CAM_INIDAT_TO_MPAS'
block => domain % blocklist
+ surface_pressure => block % time_levs(1) % state % surface_pressure % array
+ theta => 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 T = ', minval(T), maxval(T)
+ write(0,*) 'min/max for U = ', minval(U), maxval(U)
+ write(0,*) 'min/max for Tracer = ', minval(Tracer), maxval(Tracer)
+
!
- ! Fill in MPAS fields in block % time_levs(1) % state arrays from arguments
+ ! 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
-! TEMPORARY: until we are actually given fields, just call one of the test case
-! setup routines to initialize the MPAS prognostic fields
- call mpas_setup_test_case(domain)
+ !
+ ! 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(k) + Bc(k)
+ znwc(Plev+1-k+1) = Ae(k) + Be(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(Plev+1-k)
+ 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
+ theta(:,:) = 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
+ pressure(k,iCell) = pressure(k+1,iCell)-dnw(k)*h(k,iCell)
+ end do
+
+ do k=1,Plev
+ theta(k,iCell) = T(iCell,k) * (0.5*(pressure(k,iCell) + pressure(k+1,iCell))/p0)**(-rgas/cp)
+ alpha(k,iCell) = (rgas/p0) * theta(k,iCell) * (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(:), &
@@ -147,7 +301,12 @@
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
+
end subroutine cam_inidat_to_mpas
@@ -399,10 +558,36 @@
itimestep = 1
call mpas_timestep(domain, itimestep, dt)
+
+ ! Move time level 2 fields back into time level 1 for next time step
+ call shift_time_levels(domain)
end subroutine mpas_dyn_run
+!************** TEST ROUTINE -- NOT PART OF CAM-MPAS INTERFACE *************
+ subroutine mpas_write_output()
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE MPAS_WRITE_OUTPUT
+ !
+ ! Write MPAS output.nc file
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (io_output_object) :: output_obj
+
+ write(0,*) 'Called MPAS_WRITE_OUTPUT'
+
+ call output_state_init(output_obj, domain, "OUTPUT")
+ call output_state_for_domain(output_obj, domain, 1)
+ call output_state_finalize(output_obj, domain % dminfo)
+
+ end subroutine mpas_write_output
+!************** TEST ROUTINE -- NOT PART OF CAM-MPAS INTERFACE *************
+
+
subroutine mpas_final()
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SUBROUTINE MPAS_FINAL
Modified: branches/mpas_cam_coupling/src/driver_cam_interface/test_driver.F
===================================================================
--- branches/mpas_cam_coupling/src/driver_cam_interface/test_driver.F 2010-03-31 20:43:23 UTC (rev 175)
+++ branches/mpas_cam_coupling/src/driver_cam_interface/test_driver.F 2010-03-31 22:55:43 UTC (rev 176)
@@ -8,12 +8,12 @@
integer :: mpi_procID, mpi_nprocs, mpi_ierr
character (len=32) :: filename
- integer :: nCellsSolve, maxEdges, nVertLevels, nScalars
+ integer :: nCellsSolve, maxEdges, nVertLevels, nScalars, istep
real (kind=RKIND) :: Pref, Ptop
real (kind=RKIND), dimension(26) :: Ac, Bc
real (kind=RKIND), dimension(27) :: Ae, Be
real (kind=RKIND), dimension(:), pointer :: Psd, Phis
- real (kind=RKIND), dimension(:,:), pointer :: T, Ux, Uy, Omega
+ real (kind=RKIND), dimension(:,:), pointer :: T, Ux, Uy, dPdt
real (kind=RKIND), dimension(:,:), pointer :: T_tend, Ux_tend, Uy_tend
real (kind=RKIND), dimension(:,:,:), pointer :: U, Tracer
@@ -65,7 +65,7 @@
allocate(U(nCellsSolve,maxEdges,nVertLevels))
allocate(Ux(nCellsSolve,nVertLevels))
allocate(Uy(nCellsSolve,nVertLevels))
- allocate(Omega(nCellsSolve,nVertLevels))
+ allocate(dPdt(nCellsSolve,nVertLevels))
allocate(Tracer(nCellsSolve,nVertLevels,nScalars))
allocate(T_tend(nCellsSolve,nVertLevels))
allocate(Ux_tend(nCellsSolve,nVertLevels))
@@ -93,10 +93,15 @@
call mpas_init2()
- call mpas_dyn_run()
+ do istep=1,24
+ call mpas_dyn_run()
+ end do
+ ! *** Test routine to write model state -- not part of MPAS-CAM interface
+ call mpas_write_output()
+
call mpas_to_cam(nCellsSolve, nVertLevels, nScalars, &
- Psd, Phis, T, Ux, Uy, Omega, Tracer)
+ Psd, Phis, T, Ux, Uy, dPdt, Tracer)
call cam_to_mpas(nCellsSolve, nVertLevels, nScalars, &
T_tend, Ux_tend, Uy_tend, Tracer)
@@ -114,7 +119,7 @@
deallocate(U)
deallocate(Ux)
deallocate(Uy)
- deallocate(Omega)
+ deallocate(dPdt)
deallocate(Tracer)
deallocate(T_tend)
deallocate(Ux_tend)
</font>
</pre>