<p><b>duda</b> 2010-05-04 18:32:23 -0600 (Tue, 04 May 2010)</p><p>BRANCH COMMIT<br>
<br>
Update the mpas_cam_coupling branch with respect to the trunk.<br>
<br>
<br>
M namelist.input.hyd_atmos<br>
M src/core_hyd_atmos/module_test_cases.F<br>
M src/core_hyd_atmos/Registry<br>
M src/core_hyd_atmos/module_time_integration.F<br>
M src/core_sw/module_test_cases.F<br>
M src/registry/gen_inc.c<br>
M src/framework/module_dmpar.F<br>
M Makefile<br>
</p><hr noshade><pre><font color="gray">Modified: branches/mpas_cam_coupling/Makefile
===================================================================
--- branches/mpas_cam_coupling/Makefile 2010-05-04 23:21:47 UTC (rev 246)
+++ branches/mpas_cam_coupling/Makefile 2010-05-05 00:32:23 UTC (rev 247)
@@ -2,7 +2,7 @@
MODEL_FORMULATION = -DLANL_FORMULATION
EXPAND_LEVELS = -DEXPAND_LEVELS=26
-#FILE_OFFSET = -DOFFSET64BIT
+FILE_OFFSET = -DOFFSET64BIT
#########################
# Section for Zoltan TPL
@@ -75,14 +75,37 @@
"CC = mpicc" \
"SFC = gfortran" \
"SCC = gcc" \
- "FFLAGS = -O3 -m64 -ffree-line-length-none" \
+ "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3 -m64" \
"CORE = $(CORE)" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+g95:
+ ( make all \
+ "FC = mpif90" \
+ "CC = mpicc" \
+ "SFC = g95" \
+ "SCC = gcc" \
+ "FFLAGS = -O3 -ffree-line-length-huge -r8" \
+ "CFLAGS = -O3" \
+ "LDFLAGS = -O3" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+g95-serial:
+ ( make all \
+ "FC = g95" \
+ "CC = gcc" \
+ "SFC = g95" \
+ "SCC = gcc" \
+ "FFLAGS = -O3 -ffree-line-length-huge -r8" \
+ "CFLAGS = -O3" \
+ "LDFLAGS = -O3" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+
CPPINCLUDES = -I../inc -I$(NETCDF)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include
LIBS = -L$(NETCDF)/lib -lnetcdf
Modified: branches/mpas_cam_coupling/namelist.input.hyd_atmos
===================================================================
--- branches/mpas_cam_coupling/namelist.input.hyd_atmos 2010-05-04 23:21:47 UTC (rev 246)
+++ branches/mpas_cam_coupling/namelist.input.hyd_atmos 2010-05-05 00:32:23 UTC (rev 247)
@@ -13,6 +13,7 @@
config_v_theta_eddy_visc2 = 0.0
config_theta_adv_order = 2
config_scalar_adv_order = 2
+ config_mp_physics = 0
/
&io
Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/Registry
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/Registry 2010-05-04 23:21:47 UTC (rev 246)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/Registry 2010-05-05 00:32:23 UTC (rev 247)
@@ -17,6 +17,7 @@
namelist integer sw_model config_scalar_adv_order 2
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic true
+namelist integer sw_model config_mp_physics 0
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -33,10 +34,10 @@
dim maxEdges2 maxEdges2
dim nVertices nVertices
dim TWO 2
-dim R3 3
dim vertexDegree vertexDegree
dim FIFTEEN 15
dim TWENTYONE 21
+dim R3 3
dim nVertLevels nVertLevels
#dim nTracers nTracers
dim nVertLevelsP1 nVertLevels+1
Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/module_test_cases.F
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/module_test_cases.F 2010-05-04 23:21:47 UTC (rev 246)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/module_test_cases.F 2010-05-05 00:32:23 UTC (rev 247)
@@ -28,10 +28,11 @@
write(0,*) ' need hydrostatic test case configuration, error stop '
stop
- else if ((config_test_case == 1) .or. (config_test_case == 2)) then
+ else if ((config_test_case == 1) .or. (config_test_case == 2) .or. (config_test_case == 3)) then
write(0,*) ' Jablonowski and Williamson baroclinic wave test case '
if (config_test_case == 1) write(0,*) ' no initial perturbation '
if (config_test_case == 2) write(0,*) ' initial perturbation included '
+ if (config_test_case == 3) write(0,*) ' normal-mode perturbation included '
block_ptr => domain % blocklist
do while (associated(block_ptr))
call hyd_test_case_1(block_ptr % mesh, block_ptr % time_levs(1) % state, config_test_case)
@@ -69,6 +70,8 @@
real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
real (kind=RKIND), parameter :: theta_c = pii/4.0
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
+ real (kind=RKIND), parameter :: rh_max = 0.4 ! Maximum relative humidity
+ real (kind=RKIND), parameter :: k_x = 9. ! Normal mode wave number
real (kind=RKIND), dimension(:), pointer :: rdnu, rdnw, fnm, fnp, dbn, dnu, dnw
real (kind=RKIND), dimension(:), pointer :: surface_pressure
@@ -79,7 +82,12 @@
real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
real (kind=RKIND) :: ptop, p0, phi
+ real (kind=RKIND) :: lon_Edge
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature
+ real (kind=RKIND) :: ptmp, es, qvs
+ integer :: iter
+
! real (kind=RKIND), dimension(grid % nVertLevelsP1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
! real (kind=RKIND), dimension(grid % nVertLevelsP1 ) :: znuc, znuv, bn, divh, dpn, teta, phi
real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
@@ -141,13 +149,14 @@
h => state % h % array
scalars => state % scalars % array
- scalars = 0.
+ scalars(:,:,:) = 0.
p0 = 100000.
bn (1) = 1.
znw(1) = 1.
znwc(1) = 1.
- znwv(1) = (znwc(1)-.252)*pii/2.
+ !znwv(1) = (znwc(1)-.252)*pii/2.
+ znwv(1) = ((znwc(1)-.252)*pii/2.*p0-ptop)/(p0-ptop)
if (cam26) then
@@ -208,8 +217,10 @@
!
do k=1,nz1
- znuv(k ) = (znuc(k )-.252)*pii/2.
- znwv(k+1) = (znwc(k+1)-.252)*pii/2.
+ !znuv(k ) = (znuc(k )-.252)*pii/2.
+ !znwv(k+1) = (znwc(k+1)-.252)*pii/2.
+ znuv(k ) = ((znuc(k )-.252)*pii/2.*p0-ptop)/(p0-ptop)
+ znwv(k+1) = ((znwc(k+1)-.252)*pii/2.*p0-ptop)/(p0-ptop)
dnw (k) = znw(k+1)-znw(k)
rdnw(k) = 1./dnw(k)
dbn (k) = rdnw(k)*(bn(k+1)-bn(k))
@@ -247,6 +258,11 @@
r_pert = sphere_distance( grid % latEdge % array (iEdge), grid % lonEdge % array (iEdge), &
lat_pert, lon_pert, 1.)/(pert_radius)
u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1)*a/grid % dvEdge % array(iEdge)
+
+ else if (config_test_case == 3) then
+ lon_Edge = grid % lonEdge % array(iEdge)
+ u_pert = u_perturbation*cos(k_x*(lon_Edge - lon_pert)) &
+ *(0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1)))*a/grid % dvEdge % array(iEdge)
else
u_pert = 0.0
end if
@@ -274,14 +290,16 @@
)
end do
- do k=1,nz1
- if(znuc(k).ge.eta_t) then
- teta(k) = t0*znuc(k)**(rgas*dtdz/gravity)
- else
- teta(k) = t0*znuc(k)**(rgas*dtdz/gravity) + delta_t*(eta_t-znuc(k))**5
- end if
- write(6,*) ' k, reference t ',k,teta(k)
- end do
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! To get hydrostatic balance with misture -- soln. 2.
+! original scheme by Jablonowski
+! T' = -1./R_d *(p/p_0) * d(phi')/d(eta)
+! = -1./R_d * eta * d(phi')/d(eta)
+! soln. 2 -> derive temperature profile from hydrostatic balance with moisture
+!
+! T_v = -1/(1+q_v)*(p/R_d)* d(eta)/d(p_d) * d(phi)/d(eta)
+! phi'(k) = phi(k+1) + d(eta)* alpha_pert * d(eta)/d(p_d)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do iCell=1,grid % nCells
@@ -299,8 +317,13 @@
end do
do k=1,nz1
-
- theta (k,iCell) = teta(k)+.75*znuc(k)*pii*u0/rgas*sin(znuv(k)) &
+ ptmp = 0.5*(pressure(k,iCell)+pressure(k+1,iCell))
+ if (znuc(k) >= eta_t) then
+ teta(k) = t0*(ptmp/p0)**(rgas*dtdz/gravity)
+ else
+ teta(k) = t0*(ptmp/p0)**(rgas*dtdz/gravity) + delta_t*(eta_t-ptmp/p0)**5
+ end if
+ theta (k,iCell) = teta(k)+.75*ptmp/h(k,iCell)*pii*u0/rgas*sin(znuv(k)) &
*sqrt(cos(znuv(k)))* &
((-2.*sin(phi)**6 &
*(cos(phi)**2+1./3.)+10./63.) &
@@ -333,9 +356,89 @@
end do
end do
+ write(6,*) 'ptop_dry = ',ptop,' zt_dry = ',geopotential(nz,1)/gravity
+
+ write(6,*) ' full sounding for dry'
+ do k=1,nz1
+ write(6,*) k, geopotential(k,1)/gravity, 0.01*pressure(k,1), theta(k,1), &
+ theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
+ end do
+
+!
+! initialization for moisture
+!
+ if (config_mp_physics /= 0) then
+
+ do iCell=1,grid % nCells
+ do k=1,nz1
+ ptmp = 0.5*(pressure(k,iCell) + pressure(k+1,iCell))
+ if (ptmp < 50000.) then
+ rel_hum(k,iCell) = 0.0
+ else
+ rel_hum(k,iCell) = (1.-((p0-ptmp)/50000.)**1.25)
+ end if
+ rel_hum(k,iCell) = min(rh_max,rel_hum(k,iCell))
+ end do
+ end do
+ else
+ rel_hum(:,:) = 0.
+ end if
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! iteration
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ do iter=1,30
+ do iCell=1,grid % nCells
+
+ phi = grid % latCell % array (iCell)
+ do k=1,nz1
+ ptmp = 0.5*(pressure(k+1,iCell)+pressure(k,iCell))
+
+ if(znuc(k) >= eta_t) then
+ teta(k) = t0*(ptmp/p0)**(rgas*dtdz/gravity)
+ else
+ teta(k) = t0*(ptmp/p0)**(rgas*dtdz/gravity) + delta_t*(eta_t-ptmp/p0)**5
+ end if
+
+ temperature (k,iCell) = teta(k)+.75*ptmp/h(k,iCell)*pii*u0/rgas*sin(znuv(k)) &
+ *sqrt(cos(znuv(k)))* &
+ ((-2.*sin(phi)**6 &
+ *(cos(phi)**2+1./3.)+10./63.) &
+ *2.*u0*cos(znuv(k))**1.5 &
+ +(1.6*cos(phi)**3 &
+ *(sin(phi)**2+2./3.)-pii/4.)*a*omega)
+
+ temperature(k,iCell) = temperature(k,iCell)/(1.+0.61*scalars(index_qv,k,iCell))
+
+ theta (k,iCell) = temperature(k,iCell)* &
+ (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**(-rgas/cp)
+ alpha (k,iCell) = (rgas/p0)*theta(k,iCell)*(1.+1.61*scalars(index_qv,k,iCell))* &
+ (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**cvpm
+
+ if (temperature(k,iCell) > 273.15) then
+ es = 1000.*0.6112*exp(17.67*(temperature(k,iCell)-273.15)/(temperature(k,iCell)-29.65))
+ else
+ es = 1000.*0.6112*exp(21.8745584*(temperature(k,iCell)-273.16)/(temperature(k,iCell)-7.66))
+ end if
+ qvs = (287.04/461.6)*es/(ptmp-es)
+! qvs = 380.*exp(17.27*(temperature(k,iCell)-273.)/(temperature(k,iCell)-36.))/ptmp
+
+ scalars(index_qv,k,iCell) = rel_hum(k,iCell)*qvs
+ end do
+
+ do k=nz1,1,-1
+ pressure(k,iCell) = pressure(k+1,iCell)-dnw(k)*h(k,iCell)*(1.+scalars(index_qv,k,iCell))
+ geopotential(k,iCell) = geopotential(k+1,iCell)+dnw(k)*h(k,iCell)*alpha(k,iCell)
+ end do
+
+ end do
+ end do
+
write(6,*) 'ptop = ',ptop,' zt = ',geopotential(nz,1)/gravity
- write(6,*) ' full sounding '
+ write(6,*) ' full sounding with moisture'
do k=1,nz1
write(6,*) k, geopotential(k,1)/gravity, 0.01*pressure(k,1), theta(k,1), &
theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F 2010-05-04 23:21:47 UTC (rev 246)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F 2010-05-05 00:32:23 UTC (rev 247)
@@ -119,8 +119,8 @@
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % cqu % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
Modified: branches/mpas_cam_coupling/src/core_sw/module_test_cases.F
===================================================================
--- branches/mpas_cam_coupling/src/core_sw/module_test_cases.F 2010-05-04 23:21:47 UTC (rev 246)
+++ branches/mpas_cam_coupling/src/core_sw/module_test_cases.F 2010-05-05 00:32:23 UTC (rev 247)
@@ -508,7 +508,7 @@
real (kind=RKIND), intent(in) :: theta
AA = 0.5 * w * (2.0 * omega + w) * cos(theta)**2.0 + &
- 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*R**2.0 * cos(theta)**-2.0)
+ 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*R**2.0 * cos(theta)**(-2.0))
end function AA
Modified: branches/mpas_cam_coupling/src/framework/module_dmpar.F
===================================================================
--- branches/mpas_cam_coupling/src/framework/module_dmpar.F 2010-05-04 23:21:47 UTC (rev 246)
+++ branches/mpas_cam_coupling/src/framework/module_dmpar.F 2010-05-05 00:32:23 UTC (rev 247)
@@ -121,6 +121,25 @@
end subroutine dmpar_abort
+ subroutine dmpar_global_abort(mesg)
+
+ implicit none
+
+ character (len=*), intent(in) :: mesg
+
+#ifdef _MPI
+ integer :: mpi_ierr, mpi_errcode
+
+ write(0,*) trim(mesg)
+ call MPI_Abort(MPI_COMM_WORLD, mpi_errcode, mpi_ierr)
+#endif
+
+ write(0,*) trim(mesg)
+ stop
+
+ end subroutine dmpar_global_abort
+
+
subroutine dmpar_bcast_int(dminfo, i)
implicit none
Modified: branches/mpas_cam_coupling/src/registry/gen_inc.c
===================================================================
--- branches/mpas_cam_coupling/src/registry/gen_inc.c 2010-05-04 23:21:47 UTC (rev 246)
+++ branches/mpas_cam_coupling/src/registry/gen_inc.c 2010-05-05 00:32:23 UTC (rev 247)
@@ -1548,7 +1548,7 @@
/*
- * Generate code to write 1d, 2d, 3d time-invariant fields
+ * Generate code to write 0d, 1d, 2d, 3d time-invariant fields
*/
for(j=0; j<2; j++) {
for(i=0; i<=3; i++) {
</font>
</pre>