<p><b>dwj07@fsu.edu</b> 2012-09-07 09:07:43 -0600 (Fri, 07 Sep 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Merging trunk to branch.<br>
Mostly just to pick up Makefile changes related to TAU.<br>
</p><hr noshade><pre><font color="gray">Index: branches/omp_blocks/openmp_test
===================================================================
--- branches/omp_blocks/openmp_test 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test 2012-09-07 15:07:43 UTC (rev 2145)
Property changes on: branches/omp_blocks/openmp_test
___________________________________________________________________
Modified: svn:mergeinfo
## -20,3 +20,4 ##
/branches/omp_blocks/multiple_blocks:1803-2084
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+/trunk/mpas:2107-2144
\ No newline at end of property
Modified: branches/omp_blocks/openmp_test/Makefile
===================================================================
--- branches/omp_blocks/openmp_test/Makefile 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/Makefile 2012-09-07 15:07:43 UTC (rev 2145)
@@ -278,6 +278,15 @@
OMP_MESSAGE="OpenMP is off."
endif
+ifeq "$(TAU)" "true"
+ LINKER=tau_f90.sh
+ CPPINCLUDES += -DMPAS_TAU
+ TAU_MESSAGE="TAU Hooks are on."
+else
+ LINKER=$(FC)
+ TAU_MESSAGE="TAU Hooks are off."
+endif
+
ifneq ($(wildcard $(NETCDF)/lib/libnetcdff.*), ) # CHECK FOR NETCDF4
LIBS += -lnetcdff
endif # CHECK FOR NETCDF4
@@ -302,6 +311,7 @@
CC="$(CC)" \
SFC="$(SFC)" \
SCC="$(SCC)" \
+ LINKER="$(LINKER)" \
CFLAGS="$(CFLAGS)" \
FFLAGS="$(FFLAGS)" \
LDFLAGS="$(LDFLAGS)" \
@@ -318,6 +328,7 @@
@echo $(SERIAL_MESSAGE)
@echo $(PAPI_MESSAGE)
@echo $(OMP_MESSAGE)
+ @echo $(TAU_MESSAGE)
clean:
cd src; $(MAKE) clean RM="$(RM)" CORE="$(CORE)"
$(RM) $(CORE)_model.exe
@@ -349,9 +360,10 @@
@cd src; ls -d core_* | grep ".*" | sed "s/core_/ /g"
@echo ""
@echo "Available Options:"
- @echo " SERIAL=true - builds serial version. Default is parallel version."
- @echo " DEBUG=true - builds debug version. Default is optimized version."
- @echo " USE_PAPI=true - builds version using PAPI for timers and hardware counters. Default is off."
+ @echo " SERIAL=true - builds serial version. Default is parallel version."
+ @echo " DEBUG=true - builds debug version. Default is optimized version."
+ @echo " USE_PAPI=true - builds version using PAPI for timers. Default is off."
+ @echo " TAU=true - builds version using TAU hooks for profiling. Default is off."
@echo ""
@echo "Ensure that NETCDF (and PAPI if USE_PAPI=true) are environment variables"
@echo "that point to the absolute paths for the libraries."
Modified: branches/omp_blocks/openmp_test/namelist.input.init_nhyd_atmos
===================================================================
--- branches/omp_blocks/openmp_test/namelist.input.init_nhyd_atmos 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/namelist.input.init_nhyd_atmos 2012-09-07 15:07:43 UTC (rev 2145)
@@ -5,6 +5,12 @@
config_stop_time = '2010-10-23_00:00:00'
/
+&dcmip
+ config_dcmip_case = '2-0-0'
+ config_planet_scale = 1.0
+ config_rotation_rate_scale = 1.0
+/
+
&dimensions
config_nvertlevels = 41
config_nsoillevels = 4
@@ -33,8 +39,8 @@
/
&io
- config_input_name = 'x1.40962.geogrid.nc'
- config_output_name = 'x1.40962.init.2010-10-23.nc'
+ config_input_name = 'x1.40962.grid.nc'
+ config_output_name = 'x1.40962.init.nc'
config_pio_num_iotasks = 0
config_pio_stride = 1
/
@@ -47,7 +53,4 @@
/
&restart
- config_restart_interval = 3000
- config_do_restart = .false.
- config_restart_time = 1036800.0
/
Modified: branches/omp_blocks/openmp_test/namelist.input.nhyd_atmos
===================================================================
--- branches/omp_blocks/openmp_test/namelist.input.nhyd_atmos 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/namelist.input.nhyd_atmos 2012-09-07 15:07:43 UTC (rev 2145)
@@ -4,12 +4,12 @@
config_start_time = '2010-10-23_00:00:00'
config_run_duration = '5_00:00:00'
config_number_of_sub_steps = 6
- config_h_mom_eddy_visc2 = 0000.
- config_h_mom_eddy_visc4 = 0.
- config_v_mom_eddy_visc2 = 00.0
- config_h_theta_eddy_visc2 = 0000.
- config_h_theta_eddy_visc4 = 00.
- config_v_theta_eddy_visc2 = 00.0
+ config_h_mom_eddy_visc2 = 0.0
+ config_h_mom_eddy_visc4 = 0.0
+ config_v_mom_eddy_visc2 = 0.0
+ config_h_theta_eddy_visc2 = 0.0
+ config_h_theta_eddy_visc4 = 0.0
+ config_v_theta_eddy_visc2 = 0.0
config_horiz_mixing = '2d_smagorinsky'
config_len_disp = 120000.0
config_theta_adv_order = 3
@@ -35,12 +35,8 @@
config_xnutr = 0.0
/
-&dimensions
- config_nvertlevels = 41
-/
-
&io
- config_input_name = 'x1.40962.init.2010-10-23.nc'
+ config_input_name = 'x1.40962.init.nc'
config_output_name = 'x1.40962.output.nc'
config_restart_name = 'restart.nc'
config_output_interval = '1_00:00:00'
Modified: branches/omp_blocks/openmp_test/namelist.input.nhyd_atmos_jw
===================================================================
--- branches/omp_blocks/openmp_test/namelist.input.nhyd_atmos_jw 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/namelist.input.nhyd_atmos_jw 2012-09-07 15:07:43 UTC (rev 2145)
@@ -1,18 +1,17 @@
&nhyd_model
- config_test_case = 2
config_time_integration = 'SRK3'
config_dt = 450
- config_ntimesteps = 1920
- config_output_interval = 192
+ config_start_time = '0000-01-01_00:00:00'
+ config_run_duration = '10_00:00:00'
config_number_of_sub_steps = 6
- config_h_mom_eddy_visc2 = 0.0e+04
- config_h_mom_eddy_visc4 = 0.
- config_v_mom_eddy_visc2 = 00.0
- config_h_theta_eddy_visc2 = 0.0e+04
- config_h_theta_eddy_visc4 = 00.
- config_v_theta_eddy_visc2 = 00.0
+ config_h_mom_eddy_visc2 = 0.0
+ config_h_mom_eddy_visc4 = 0.0
+ config_v_mom_eddy_visc2 = 0.0
+ config_h_theta_eddy_visc2 = 0.0
+ config_h_theta_eddy_visc4 = 0.0
+ config_v_theta_eddy_visc2 = 0.0
config_horiz_mixing = '2d_smagorinsky'
- config_len_disp = 60000.
+ config_len_disp = 120000.
config_u_vadv_order = 3
config_w_vadv_order = 3
config_theta_vadv_order = 3
@@ -39,7 +38,7 @@
/
&io
- config_input_name = 'grid.nc'
+ config_input_name = 'x1.40962.init.nc'
config_output_name = 'output.nc'
config_restart_name = 'restart.nc'
config_pio_num_iotasks = 0
@@ -48,15 +47,14 @@
&decomposition
config_number_of_blocks = 0
- config_block_decomp_file_prefix = 'graph.info.part.'
+ config_block_decomp_file_prefix = 'x1.40962.graph.info.part.'
config_explicit_proc_decomp = .false.
config_proc_decomp_file_prefix = 'graph.info.part.'
/
&restart
- config_restart_interval = 3000
+ config_restart_interval = '10_00:00:00'
config_do_restart = .false.
- config_restart_time = 1036800.0
/
&physics
Modified: branches/omp_blocks/openmp_test/src/Makefile
===================================================================
--- branches/omp_blocks/openmp_test/src/Makefile 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/Makefile 2012-09-07 15:07:43 UTC (rev 2145)
@@ -3,7 +3,7 @@
all: mpas
mpas: reg_includes externals frame ops dycore drver
- $(FC) $(LDFLAGS) -o $(CORE)_model.exe driver/*.o -L. -ldycore -lops -lframework $(LIBS) -I./external/esmf_time_f90 -L./external/esmf_time_f90 -lesmf_time
+ $(LINKER) $(LDFLAGS) -o $(CORE)_model.exe driver/*.o -L. -ldycore -lops -lframework $(LIBS) -I./external/esmf_time_f90 -L./external/esmf_time_f90 -lesmf_time
reg_includes:
( cd registry; $(MAKE) CC="$(SCC)" )
Modified: branches/omp_blocks/openmp_test/src/core_hyd_atmos/mpas_atmh_time_integration.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_hyd_atmos/mpas_atmh_time_integration.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_hyd_atmos/mpas_atmh_time_integration.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -276,7 +276,7 @@
block % mesh % areaCell % array (iCell) &
- block % state % time_levs(2) % state % pressure % array (block % mesh % nVertLevels + 1, 1) * &
block % mesh % areaCell % array (iCell)
- do k=1, block % mesh % nVertLevelsSolve
+ do k=1, block % mesh % nVertLevels ! Could be nVertLevelsSolve?
scalar_mass = scalar_mass - block % state % time_levs(2) % state % scalars % array (2,k,iCell) * &
block % state % time_levs(2) % state % h % array (k,iCell) * &
block % mesh % dnw % array (k) * &
@@ -1378,7 +1378,7 @@
end do
wdtn(:,nVertLevels+1) = 0.
- do k=1,grid % nVertLevelsSolve
+ do k=1,grid % nVertLevels ! Could be nVertLevelsSolve?
do iScalar=1,num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
Modified: branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/Registry 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/Registry 2012-09-07 15:07:43 UTC (rev 2145)
@@ -8,6 +8,9 @@
namelist integer nhyd_model config_theta_adv_order 3
namelist real nhyd_model config_coef_3rd_order 0.25
namelist integer nhyd_model config_num_halos 2
+namelist character dcmip config_dcmip_case 2-0-0
+namelist real dcmip config_planet_scale 1.0
+namelist real dcmip config_rotation_rate_scale 1.0
namelist integer dimensions config_nvertlevels 26
namelist integer dimensions config_nsoillevels 4
namelist integer dimensions config_nfglevels 27
Modified: branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/mpas_init_atm_mpas_core.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/mpas_init_atm_mpas_core.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/mpas_init_atm_mpas_core.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -22,6 +22,7 @@
block => domain % blocklist
do while (associated(block))
block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
+ block % mesh % sphere_radius = a / config_planet_scale
block => block % next
end do
Modified: branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -121,14 +121,61 @@
block_ptr => block_ptr % next
end do
+ else if (config_test_case == 9 ) then
+
+ write(0,*) ' '
+ write(0,*) ' '
+ write(0,*) ' Setting up DCMIP test case '//trim(config_dcmip_case)
+ write(0,*) ' '
+ write(0,*) ' '
+
+ if (trim(config_dcmip_case) == '2-0-0' .or. &
+ trim(config_dcmip_case) == '2-0-1') then
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call init_atm_test_case_resting_atmosphere(block_ptr % mesh, block_ptr % state % time_levs(1) % state, &
+ block_ptr % diag, config_test_case)
+ block_ptr => block_ptr % next
+ end do
+
+ else if (trim(config_dcmip_case) == '2-1' .or. &
+ trim(config_dcmip_case) == '2-1a' .or. &
+ trim(config_dcmip_case) == '2-2' .or. &
+ trim(config_dcmip_case) == '3-1') then
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call init_atm_test_case_reduced_radius(block_ptr % mesh, block_ptr % state % time_levs(1) % state, &
+ block_ptr % diag, config_test_case)
+ block_ptr => block_ptr % next
+ end do
+
+ else
+
+ write(0,*) ' '
+ write(0,*) ' *************'
+ write(0,*) ' Unrecognized DCMIP case '//trim(config_dcmip_case)
+ write(0,*) ' Please choose either 2-0-0, 2-0-1, 2-1, 2-1a, 2-2, or 3-1'
+ write(0,*) ' *************'
+ write(0,*) ' '
+ call mpas_dmpar_abort(domain % dminfo)
+
+ end if
+
else
- write(0,*) ' Only test cases 1, 2, 3, 4, 5, 6, 7, and 8 are currently supported for nonhydrostatic core '
- stop
+ write(0,*) ' '
+ write(0,*) ' *************'
+ write(0,*) ' Only test cases 1 through 9 are currently supported for the nonhydrostatic core'
+ write(0,*) ' *************'
+ write(0,*) ' '
+ call mpas_dmpar_abort(domain % dminfo)
end if
+ ! Copy initialized state to all time levels
block_ptr => domain % blocklist
do while (associated(block_ptr))
do i=2,nTimeLevs
@@ -166,7 +213,10 @@
real (kind=RKIND), parameter :: u0 = 35.0
real (kind=RKIND), parameter :: alpha_grid = 0. ! no grid rotation
- real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+
+! real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+ real (kind=RKIND) :: omega_e
+
real (kind=RKIND), parameter :: t0b = 250., t0 = 288., delta_t = 4.8e+05, dtdz = 0.005, eta_t = 0.2
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
@@ -229,23 +279,25 @@
real (kind=RKIND) :: z_edge, z_edge3, d2fdx2_cell1, d2fdx2_cell2
logical, parameter :: moisture = .true.
+! logical, parameter :: moisture = .false.
+
!
- ! Scale all distances and areas from a unit sphere to one with radius a
+ ! Scale all distances and areas from a unit sphere to one with radius sphere_radius
!
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
+ grid % xCell % array = grid % xCell % array * grid % sphere_radius
+ grid % yCell % array = grid % yCell % array * grid % sphere_radius
+ grid % zCell % array = grid % zCell % array * grid % sphere_radius
+ grid % xVertex % array = grid % xVertex % array * grid % sphere_radius
+ grid % yVertex % array = grid % yVertex % array * grid % sphere_radius
+ grid % zVertex % array = grid % zVertex % array * grid % sphere_radius
+ grid % xEdge % array = grid % xEdge % array * grid % sphere_radius
+ grid % yEdge % array = grid % yEdge % array * grid % sphere_radius
+ grid % zEdge % array = grid % zEdge % array * grid % sphere_radius
+ grid % dvEdge % array = grid % dvEdge % array * grid % sphere_radius
+ grid % dcEdge % array = grid % dcEdge % array * grid % sphere_radius
+ grid % areaCell % array = grid % areaCell % array * grid % sphere_radius**2.0
+ grid % areaTriangle % array = grid % areaTriangle % array * grid % sphere_radius**2.0
+ grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * grid % sphere_radius**2.0
weightsOnEdge => grid % weightsOnEdge % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
@@ -311,7 +363,8 @@
znut = eta_t
etavs = (1.-0.252)*pii/2.
- r_earth = a
+ r_earth = grid % sphere_radius
+ omega_e = omega * config_rotation_rate_scale
p0 = 1.e+05
write(0,*) ' point 1 in test case setup '
@@ -518,9 +571,14 @@
ppi(1) = ppi(1)-ppb(1,i)
do k=1,nz1-1
- ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity* &
- (rr(k ,i)+(rr(k ,i)+rb(k ,i))*qv_2d(k ,i) &
- +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*qv_2d(k+1,i))
+
+! ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity* &
+! (rr(k ,i)+(rr(k ,i)+rb(k ,i))*qv_2d(k ,i) &
+! +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*qv_2d(k+1,i))
+
+ ppi(k+1) = ppi(k)-dzu(k+1)*gravity* &
+ ( (rr(k ,i)+(rr(k ,i)+rb(k ,i))*qv_2d(k ,i))*fzp(k+1) &
+ + (rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*qv_2d(k+1,i))*fzm(k+1) )
end do
do k=1,nz1
@@ -550,7 +608,7 @@
end do
call init_atm_recompute_geostrophic_wind(u_2d,rho_2d,pp_2d,qv_2d,lat_2d,zz_2d,zx_2d, &
- cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat)
+ cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat,grid%sphere_radius)
end if
@@ -651,9 +709,15 @@
ppi(1) = ppi(1)-ppb(1,i)
do k=1,nz1-1
- ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity* &
- (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i) &
- +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))
+
+! ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity* &
+! (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i) &
+! +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))
+
+ ppi(k+1) = ppi(k)-dzu(k+1)*gravity* &
+ ( (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i))*fzp(k+1) &
+ + (rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))*fzm(k+1) )
+
end do
do k=1,nz1
@@ -701,24 +765,24 @@
lat2 = grid%latVertex%array(vtx2)
iCell1 = grid % cellsOnEdge % array(1,iEdge)
iCell2 = grid % cellsOnEdge % array(2,iEdge)
- flux = (0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1)))*a/grid % dvEdge % array(iEdge)
+ flux = (0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1))) * grid % sphere_radius / grid % dvEdge % array(iEdge)
if (config_test_case == 2) then
r_pert = sphere_distance( grid % latEdge % array (iEdge), grid % lonEdge % array (iEdge), &
lat_pert, lon_pert, 1.0_RKIND)/(pert_radius)
- u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1)*a/grid % dvEdge % array(iEdge)
+ u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1) * grid % sphere_radius / 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)
+ *(0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1))) * grid % sphere_radius / grid % dvEdge % array(iEdge)
else
u_pert = 0.0
end if
if (rebalance) then
- call init_atm_calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1,lat2,grid % dvEdge % array(iEdge),a,u0,nz1,nlat)
+ call init_atm_calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1,lat2,grid % dvEdge % array(iEdge),grid%sphere_radius,u0,nz1,nlat)
do k=1,grid % nVertLevels
fluxk = u0*flux_zonal(k)/(0.5*(rb(k,iCell1)+rb(k,iCell2)+rr(k,iCell1)+rr(k,iCell2)))
state % u % array(k,iEdge) = fluxk + u_pert
@@ -744,14 +808,14 @@
! Generate rotated Coriolis field
!
- grid % fEdge % array(iEdge) = 2.0 * omega * &
+ grid % fEdge % array(iEdge) = 2.0 * omega_e * &
( -cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha_grid) + &
sin(grid%latEdge%array(iEdge)) * cos(alpha_grid) &
)
end do
do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 2.0 * omega * &
+ grid % fVertex % array(iVtx) = 2.0 * omega_e * &
(-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha_grid) + &
sin(grid%latVertex%array(iVtx)) * cos(alpha_grid) &
)
@@ -778,11 +842,18 @@
d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
- d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
- end do
+! WCS fix 20120711
+
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
- (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
@@ -871,6 +942,7 @@
end subroutine init_atm_test_case_jw
+
subroutine init_atm_calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1_in,lat2_in,dvEdge,a,u0,nz1,nlat)
implicit none
@@ -920,9 +992,11 @@
end subroutine init_atm_calc_flux_zonal
+
+
!SHP-balance
subroutine init_atm_recompute_geostrophic_wind(u_2d,rho_2d,pp_2d,qv_2d,lat_2d,zz_2d,zx_2d, &
- cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat)
+ cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat,rad)
implicit none
integer, intent(in) :: nz1,nlat
@@ -931,18 +1005,21 @@
real (kind=RKIND), dimension(nz1,nlat-1), intent(in) :: zx_2d
real (kind=RKIND), dimension(nlat), intent(in) :: lat_2d
real (kind=RKIND), dimension(nz1), intent(in) :: fzm, fzp, rdzw
- real (kind=RKIND), intent(in) :: cf1, cf2, cf3, dlat
+ real (kind=RKIND), intent(in) :: cf1, cf2, cf3, dlat, rad
!local variable
real (kind=RKIND), dimension(nz1,nlat-1) :: pgrad, ru, u
real (kind=RKIND), dimension(nlat-1) :: f
real (kind=RKIND), dimension(nz1+1) :: dpzx
- real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+! real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+ real (kind=RKIND) :: omega_e
+
real (kind=RKIND) :: rdx, qtot, r_earth, phi
integer :: k,i, itr
- r_earth = a
+ r_earth = rad
+ omega_e = omega * config_rotation_rate_scale
rdx = 1./(dlat*r_earth)
do i=1,nlat-1
@@ -1012,10 +1089,9 @@
u_2d(k,i) = (3.*u_2d(k,i-1)-u_2d(k,i-2))*.5
end do
-
end subroutine init_atm_recompute_geostrophic_wind
-!----------------------------------------------------------------------------------------------------------
+
subroutine init_atm_test_case_squall_line(dminfo, grid, state, diag, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup squall line and supercell test case
@@ -1999,7 +2075,7 @@
write(0,*) ' *** sounding for the simulation ***'
write(0,*) ' z theta pres qv rho_m u rr'
do k=1,nz1
- write(6,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
+ write(0,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
.01*p0*p(k,1)**(1./rcp), &
1000.*scalars(index_qv,k,1), &
@@ -2179,7 +2255,10 @@
real (kind=RKIND), parameter :: u0 = 35.0
real (kind=RKIND), parameter :: alpha_grid = 0. ! no grid rotation
- real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+
+! real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+ real (kind=RKIND) :: omega_e
+
real (kind=RKIND), parameter :: t0b = 250., t0 = 288., delta_t = 4.8e+05, dtdz = 0.005, eta_t = 0.2
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
@@ -2338,7 +2417,8 @@
etavs = (1.-0.252)*pii/2.
rcv = rgas/(cp-rgas)
- r_earth = a
+ r_earth = grid % sphere_radius
+ omega_e = omega
p0 = 1.e+05
interp_list(1) = FOUR_POINT
@@ -2347,25 +2427,25 @@
!
- ! Scale all distances and areas from a unit sphere to one with radius a
+ ! Scale all distances and areas from a unit sphere to one with radius sphere_radius
!
if (config_static_interp) then
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
+ grid % xCell % array = grid % xCell % array * r_earth
+ grid % yCell % array = grid % yCell % array * r_earth
+ grid % zCell % array = grid % zCell % array * r_earth
+ grid % xVertex % array = grid % xVertex % array * r_earth
+ grid % yVertex % array = grid % yVertex % array * r_earth
+ grid % zVertex % array = grid % zVertex % array * r_earth
+ grid % xEdge % array = grid % xEdge % array * r_earth
+ grid % yEdge % array = grid % yEdge % array * r_earth
+ grid % zEdge % array = grid % zEdge % array * r_earth
+ grid % dvEdge % array = grid % dvEdge % array * r_earth
+ grid % dcEdge % array = grid % dcEdge % array * r_earth
+ grid % areaCell % array = grid % areaCell % array * r_earth**2.0
+ grid % areaTriangle % array = grid % areaTriangle % array * r_earth**2.0
+ grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * r_earth**2.0
scalars(:,:,:) = 0.
@@ -3204,9 +3284,7 @@
tempField % prev => null()
tempField % next => null()
- call mpas_timer_start("EXCHANGE_1D_REAL")
call mpas_dmpar_exch_halo_field(tempField)
- call mpas_timer_stop("EXCHANGE_1D_REAL")
! dzmina = minval(hs(:)-hx(k-1,:))
dzmina = minval(zw(k)+ah(k)*hs(1:grid%nCellsSolve)-zw(k-1)-ah(k-1)*hx(k-1,1:grid%nCellsSolve))
@@ -4308,6 +4386,7 @@
end subroutine init_atm_test_case_gfs
+
subroutine init_atm_test_case_sfc(domain, dminfo, grid, fg, state, diag, test_case, parinfo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Real-data test case using SST data
@@ -4544,6 +4623,1583 @@
end subroutine init_atm_test_case_sfc
+!--------------------- TEST CASE 9 -----------------------------------------------
+
+
+ subroutine init_atm_test_case_reduced_radius(grid, state, diag, test_case)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Setup Schar-type mountain wave test case on reduced radius sphere
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+ type (state_type), intent(inout) :: state
+ type (diag_type), intent(inout) :: diag
+ integer, intent(in) :: test_case
+
+ real (kind=RKIND), parameter :: t0=300., hm=250., alpha=0.
+! real (kind=RKIND), parameter :: t0=288., hm=0., alpha=0.
+
+ ! Parameters for test case 3-1
+ real (kind=RKIND), parameter :: widthParm = 5000.0, &
+ dTheta = 1.0, &
+ L_z = 20000.0, &
+ theta_c = 0.0, &
+ lambda_c = 2.0 * pii / 3.0
+
+
+ real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
+ real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
+ real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars, deriv_two, zb, zb3
+
+ !This is temporary variable here. It just need when calculate tangential velocity v.
+ integer :: eoe, j
+ integer, dimension(:), pointer :: nEdgesOnEdge, nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge, edgesOnCell
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell, dcEdge
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
+
+ integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, kz, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
+ integer :: index_qv
+
+ real (kind=RKIND) :: ptop, p0, pis, flux, d2fdx2_cell1, d2fdx2_cell2
+
+ real (kind=RKIND) :: ztemp, zd, zt, dz, str
+ real(kind=RKIND), dimension(:), pointer :: hs, hs1
+
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh
+ real (kind=RKIND) :: es, qvs, xnutr, ptemp
+ integer :: iter, nsm
+ integer, dimension(:,:), pointer :: cellsOnCell
+
+ type (field1DReal), pointer :: tempField
+ type (field1DReal), target :: tempFieldTarget
+
+ type (block_type), pointer :: block
+ type (parallel_info), pointer :: parinfo
+ type (dm_info), pointer :: dminfo
+
+ real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
+ real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex
+
+ real (kind=RKIND) :: d1, d2, d3, cof1, cof2, cf1, cf2, cf3
+ real (kind=RKIND) :: um, us, rcp, rcv
+ real (kind=RKIND) :: xmid, temp, pres, a_scale, xac, xlac, shear, tsurf, usurf
+
+ real (kind=RKIND) :: xi, yi, ri, xa, xc, yc, xla, zinv, xn2, xn2m, xn2l, sm, dzh, dzht, dzmin, dzmina, dzminf, &
+ dzmina_global, z_edge, z_edge3, sm0
+ real (kind=RKIND) :: theta_pert, s
+
+ integer, dimension(grid % nCells, 2) :: next_cell
+ real (kind=RKIND), dimension(grid % nCells) :: hxzt, pitop, ptopb
+ logical, parameter :: terrain_smooth = .false.
+
+ block => grid % block
+ parinfo => block % parinfo
+ dminfo => block % domain % dminfo
+
+
+ !
+ ! Scale all distances
+ !
+ a_scale = grid % sphere_radius
+
+ grid % xCell % array = grid % xCell % array * a_scale
+ grid % yCell % array = grid % yCell % array * a_scale
+ grid % zCell % array = grid % zCell % array * a_scale
+ grid % xVertex % array = grid % xVertex % array * a_scale
+ grid % yVertex % array = grid % yVertex % array * a_scale
+ grid % zVertex % array = grid % zVertex % array * a_scale
+ grid % xEdge % array = grid % xEdge % array * a_scale
+ grid % yEdge % array = grid % yEdge % array * a_scale
+ grid % zEdge % array = grid % zEdge % array * a_scale
+ grid % dvEdge % array = grid % dvEdge % array * a_scale
+ grid % dcEdge % array = grid % dcEdge % array * a_scale
+ grid % areaCell % array = grid % areaCell % array * a_scale**2.0
+ grid % areaTriangle % array = grid % areaTriangle % array * a_scale**2.0
+ grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a_scale**2.0
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ AreaCell => grid % AreaCell % array
+ CellsOnEdge => grid % CellsOnEdge % array
+ cellsOnCell => grid % cellsOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ deriv_two => grid % deriv_two % array
+
+ nz1 = grid % nVertLevels
+ nz = nz1 + 1
+ nCellsSolve = grid % nCellsSolve
+
+ zgrid => grid % zgrid % array
+ zb => grid % zb % array
+ zb3 => grid % zb3 % array
+ rdzw => grid % rdzw % array
+ dzu => grid % dzu % array
+ rdzu => grid % rdzu % array
+ fzm => grid % fzm % array
+ fzp => grid % fzp % array
+ zx => grid % zx % array
+ zz => grid % zz % array
+ hx => grid % hx % array
+ dss => grid % dss % array
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+
+ ppb => diag % pressure_base % array
+ pb => diag % exner_base % array
+ rb => diag % rho_base % array
+ tb => diag % theta_base % array
+ rtb => diag % rtheta_base % array
+ p => diag % exner % array
+ cqw => diag % cqw % array
+
+ rho_zz => state % rho_zz % array
+
+ pp => diag % pressure_p % array
+ rr => diag % rho_p % array
+ t => state % theta_m % array
+ rt => diag % rtheta_p % array
+ u => state % u % array
+ ru => diag % ru % array
+
+ scalars => state % scalars % array
+
+ index_qv = state % index_qv
+
+ scalars(:,:,:) = 0.
+
+ call atm_initialize_advection_rk(grid)
+ call atm_initialize_deformation_weights(grid)
+
+ if (trim(config_dcmip_case) == '2-1') then
+ zt = 30000.
+ xnutr = 0.1 ! Coefficient for implicit w damping in absorbing layer
+ zd = 20000. ! Bottom of absorbing layer
+ write(0,*) ' test case 2-1, zt, zd, xnutr ', zt,zd,xnutr
+ end if
+
+ if (trim(config_dcmip_case) == '2-1a') then
+ zt = 20000.
+ xnutr = 0.1 ! Coefficient for implicit w damping in absorbing layer
+ zd = 10000. ! Bottom of absorbing layer
+ write(0,*) ' test case 2-1a, zt, zd, xnutr ', zt,zd,xnutr
+ end if
+
+ if (trim(config_dcmip_case) == '2-2') then
+ zt = 30000.
+ xnutr = 0.1 ! Coefficient for implicit w damping in absorbing layer
+ zd = 20000. ! Bottom of absorbing layer
+ write(0,*) ' test case 2-2, zt, zd, xnutr ', zt,zd,xnutr
+ end if
+
+ if (trim(config_dcmip_case) == '3-1') then
+ zt = 10000.
+ xnutr = 0.0 ! Coefficient for implicit w damping in absorbing layer
+ zd = 10000. ! Bottom of absorbing layer
+ write(0,*) ' test case 3-1, zt, zd, xnutr ', zt,zd,xnutr
+ end if
+
+ p0 = 1.e+05
+ rcp = rgas/cp
+ rcv = rgas/(cp-rgas)
+
+ ! metrics for hybrid coordinate and vertical stretching
+ str = 1.0
+
+
+ dz = zt/float(nz1)
+! write(0,*) ' dz = ',dz
+
+ do k=1,nz
+
+! sh(k) is the stretching specified for height surfaces
+
+ zc(k) = zt*(real(k-1)*dz/zt)**str
+
+! to specify specific heights zc(k) for coordinate surfaces,
+! input zc(k)
+! zw(k) is the hieght of zeta surfaces
+! zw(k) = (k-1)*dz yields constant dzeta
+! and nonconstant dzeta/dz
+! zw(k) = sh(k)*zt yields nonconstant dzeta
+! and nearly constant dzeta/dz
+
+! zw(k) = float(k-1)*dz
+ zw(k) = zc(k)
+!
+! ah(k) governs the transition between terrain-following
+! and pureheight coordinates
+! ah(k) = 0 is a terrain-following coordinate
+! ah(k) = 1 is a height coordinate
+
+! ah(k) = 1.-cos(.5*pii*(k-1)*dz/zt)**6
+ ah(k) = 1.
+! write(0,*) ' k, zc, zw, ah ',k,zc(k),zw(k),ah(k)
+ end do
+ do k=1,nz1
+ dzw (k) = zw(k+1)-zw(k)
+ rdzw(k) = 1./dzw(k)
+ zu(k ) = .5*(zw(k)+zw(k+1))
+ end do
+ do k=2,nz1
+ dzu (k) = .5*(dzw(k)+dzw(k-1))
+ rdzu(k) = 1./dzu(k)
+ fzp (k) = .5* dzw(k )/dzu(k)
+ fzm (k) = .5* dzw(k-1)/dzu(k)
+ rdzwp(k) = dzw(k-1)/(dzw(k )*(dzw(k)+dzw(k-1)))
+ rdzwm(k) = dzw(k )/(dzw(k-1)*(dzw(k)+dzw(k-1)))
+ end do
+
+!********** how are we storing cf1, cf2 and cf3?
+
+ d1 = .5*dzw(1)
+ d2 = dzw(1)+.5*dzw(2)
+ d3 = dzw(1)+dzw(2)+.5*dzw(3)
+ !cf1 = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ !cf2 = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ !cf3 = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+
+ cof1 = (2.*dzu(2)+dzu(3))/(dzu(2)+dzu(3))*dzw(1)/dzu(2)
+ cof2 = dzu(2) /(dzu(2)+dzu(3))*dzw(1)/dzu(3)
+ cf1 = fzp(2) + cof1
+ cf2 = fzm(2) - cof1 - cof2
+ cf3 = cof2
+
+ grid % cf1 % scalar = cf1
+ grid % cf2 % scalar = cf2
+ grid % cf3 % scalar = cf3
+
+ write(0,*) 'EARTH RADIUS = ', grid % sphere_radius
+
+! setting for terrain
+
+! MGD for both 2-1 and 2-1a (and 2-2)
+ if (trim(config_dcmip_case) == '2-1' .or. &
+ trim(config_dcmip_case) == '2-1a' .or. &
+ trim(config_dcmip_case) == '2-2') then
+ xa = 5000.
+ xla = 4000.
+ end if
+
+ write(0,*) ' hm, xa, xla ',hm,xa,xla
+
+ hx = 0.
+
+ do iCell=1,grid % nCells
+
+ xi = grid % lonCell % array(iCell)
+ yi = grid % latCell % array(iCell)
+ xc = sphere_distance(yi, xi, yi, 0., grid % sphere_radius)
+ yc = sphere_distance(yi, xi, 0., xi, grid % sphere_radius)
+ xac = sphere_distance(yi, xa /grid % sphere_radius, yi, 0., grid % sphere_radius)
+ xlac = sphere_distance(yi, xla/grid % sphere_radius, yi, 0., grid % sphere_radius)
+
+ ri = sphere_distance(yi, xi, 0., 0., grid % sphere_radius)
+
+! MGD BEGIN 2-1
+! Circular mountain with Schar mtn cross section
+ if (trim(config_dcmip_case) == '2-1') then
+ hx(1,iCell) = hm*exp(-(ri/xa)**2)*cos(pii*ri/xla)**2
+ end if
+! MGD END 2-1
+
+! MGD BEGIN 2-2
+! Circular mountain with Schar mtn cross section
+ if (trim(config_dcmip_case) == '2-2') then
+ hx(1,iCell) = hm*exp(-(ri/xa)**2)*cos(pii*ri/xla)**2
+ end if
+! MGD END 2-2
+
+! MGD BEGIN 2-1a
+! proposed to be run with x333 rather than x500
+! Ridge mountain with Schar mtn cross section
+ if (trim(config_dcmip_case) == '2-1a') then
+ hx(1,iCell) = hm*exp(-(xc/xac)**2)*cos(pii*xc/xlac)**2*cos(yc/grid % sphere_radius)
+ end if
+! MGD END 2-1a
+
+ hx(nz,iCell) = zt
+
+
+ enddo
+ write(0,*) ' hx computation complete '
+
+!!! MGD WE NEED TO REPLACE THIS TERRAIN SMOOTHING WITH TC9
+
+ kz = nz
+
+ if (config_smooth_surfaces) then
+
+ write(0,*) ' '
+ write(0,*) ' Smoothing vertical coordinate surfaces'
+ write(0,*) ' '
+
+ allocate(hs (grid % nCells+1))
+ allocate(hs1(grid % nCells+1))
+
+ dzmin = 0.5
+ sm0 = 0.5
+ nsm = 30
+
+ write(6,*) 'dzmin = ',dzmin,' sm0 = ',sm0,' nsm = ',nsm
+
+ do k=2,kz-1
+ hx(k,:) = hx(k-1,:)
+ dzminf = zw(k)-zw(k-1)
+
+! dzmin = max(0.5_RKIND,1.-.5*zw(k)/hm)
+
+ sm = sm0*max( min(.5*zw(k)/hm,1.0_RKIND), .05 )
+
+ do i=1,nsm
+ do iCell=1,grid % nCells
+ hs1(iCell) = 0.
+ do j = 1,nEdgesOnCell(iCell)
+
+ hs1(iCell) = hs1(iCell) + dvEdge(edgesOnCell(j,iCell)) &
+ / dcEdge(edgesOnCell(j,iCell)) &
+ * (hx(k,cellsOnCell(j,iCell))-hx(k,iCell))
+ end do
+ hs1(iCell) = hx(k,iCell) + sm*hs1(iCell)
+
+ hs(iCell) = 0.
+ ! do j = 1,nEdgesOnCell(iCell)
+ ! hs(iCell) = hs(iCell) + dvEdge(edgesOnCell(j,iCell)) &
+ ! / dcEdge(edgesOnCell(j,iCell)) &
+ ! * (hs1(cellsOnCell(j,iCell))-hs1(iCell))
+ ! end do
+ hs(iCell) = hs1(iCell) - 0.*hs(iCell)
+
+ end do
+
+ tempField => tempFieldTarget
+ tempField % block => block
+ tempField % dimSizes(1) = grid % nCells
+ tempField % sendList => parinfo % cellsToSend
+ tempField % recvList => parinfo % cellsToRecv
+ tempField % copyList => parinfo % cellsToCopy
+ tempField % array => hs
+ tempField % prev => null()
+ tempField % next => null()
+
+ call mpas_dmpar_exch_halo_field(tempField)
+
+ ! dzmina = minval(hs(:)-hx(k-1,:))
+ dzmina = minval(zw(k)+ah(k)*hs(1:grid%nCellsSolve)-zw(k-1)-ah(k-1)*hx(k-1,1:grid%nCellsSolve))
+ call mpas_dmpar_min_real(dminfo, dzmina, dzmina_global)
+ ! write(0,*) ' k,i, dzmina, dzmin, zw(k)-zw(k-1) ', k,i, dzmina, dzmin, zw(k)-zw(k-1)
+ if (dzmina_global >= dzmin*(zw(k)-zw(k-1))) then
+ hx(k,:)=hs(:)
+ dzminf = dzmina_global
+ else
+ exit
+ end if
+ end do
+ write(0,*) k,i,sm,dzminf/(zw(k)-zw(k-1)),dzmina/(zw(k)-zw(k-1))
+ end do
+
+ do k=kz,nz
+ hx(k,:) = 0.
+ end do
+
+ deallocate(hs )
+ deallocate(hs1)
+
+ else
+
+ do k=2,nz1
+ dzmina = minval(zw(k)+ah(k)*hx(k,:)-zw(k-1)-ah(k-1)*hx(k-1,:))
+ write(0,*) k,dzmina/(zw(k)-zw(k-1))
+ end do
+
+ end if
+
+
+ do iCell=1,grid % nCells
+ do k=1,nz
+ if (config_smooth_surfaces) then
+ zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(k,iCell)/zt)+hx(k,iCell)) &
+ + (1.-ah(k)) * zc(k)
+ else
+ zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(1,iCell)/zt)+hx(1,iCell)) &
+ + (1.-ah(k)) * zc(k)
+ end if
+ end do
+ do k=1,nz1
+ zz (k,iCell) = (zw(k+1)-zw(k))/(zgrid(k+1,iCell)-zgrid(k,iCell))
+ end do
+ end do
+
+ do i=1, grid % nEdges
+ iCell1 = grid % CellsOnEdge % array(1,i)
+ iCell2 = grid % CellsOnEdge % array(2,i)
+ do k=1,nz
+ zx (k,i) = (zgrid(k,iCell2)-zgrid(k,iCell1)) / grid % dcEdge % array(i)
+ end do
+ end do
+ do i=1, grid % nCells
+ do k=1,nz1
+ ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
+ dss(k,i) = 0.
+ ztemp = zgrid(k,i)
+ if(ztemp.gt.zd+.1) then
+ dss(k,i) = dss(k,i)+xnutr*sin(.5*pii*(ztemp-zd)/(zt-zd))**2
+ end if
+ end do
+ enddo
+
+ write(0,*) ' grid metrics setup complete '
+
+!
+! mountain wave initialization
+!
+!MGD BEGIN 3-1
+! Coefficients used to initialize 2 layer sounding based on stability
+ if (trim(config_dcmip_case) == '3-1') then
+ zinv = 3000. ! Height of lower layer
+ xn2 = 0.0001 ! N^2 for upper layer
+ xn2m = 0.0001 ! N^2 for reference sounding
+ xn2l = 0.0001 ! N^@ for lower layer
+ end if
+!MGD END 3-1
+
+ if (trim(config_dcmip_case) == '2-1' .or. &
+ trim(config_dcmip_case) == '2-1a' .or. &
+ trim(config_dcmip_case) == '2-2' .or. &
+ trim(config_dcmip_case) == '3-1') then
+ um = 20. ! base wind for 2-1, 2-1a, 2-2, and 3-1
+ end if
+
+ if (trim(config_dcmip_case) == '2-2') then
+ shear = 0.00025 ! MGD 2-2
+ else
+ shear = 0. ! MGD everything else, 2-1, ...
+ end if
+
+ do i=1,grid % nCells
+
+! Surface temp and Exner function as function of latitude to balance wind fed
+
+ tsurf = t0*exp(-shear*um**2/gravity*sin(grid%latCell%array(i))**2)
+ pis = exp(-um**2*sin(grid%latCell%array(i))**2/(2.*cp*tsurf))
+
+ do k=1,nz1
+ ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
+
+!MGD FOR 2-1, 2-1a, 2-2
+! Isothermal temerature initialization
+ if (trim(config_dcmip_case) == '2-1' .or. &
+ trim(config_dcmip_case) == '2-1a' .or. &
+ trim(config_dcmip_case) == '2-2') then
+
+ t (k,i) = tsurf/pis*exp(gravity*ztemp/(cp*tsurf))
+ tb (k,i) = t0*exp(gravity*ztemp/(cp*t0))
+!! JBK fix, 20120801
+ !! tb(k,i) = t(k,i)
+
+ end if
+
+!MGD FOR 3-1
+! Initialization based on stability
+ if (trim(config_dcmip_case) == '3-1') then
+ if(ztemp .le. zinv) then
+ t (k,i) = t0*(1.+xn2l/gravity*ztemp)
+ else
+ t (k,i) = t0*(1.+xn2l/gravity*zinv+xn2/gravity*(ztemp-zinv))
+ end if
+ tb(k,i) = t0*(1. + xn2m/gravity*ztemp)
+ end if
+
+ rh(k,i) = 0.
+ end do
+
+
+! MGD ADD CODE HERE FOR 3-1 THERMAL PERT
+ if (trim(config_dcmip_case) == '3-1') then
+ do k=1,nz1
+ s = widthParm**2.0 / (widthParm**2.0 + sphere_distance(theta_c, lambda_c, &
+ grid%latCell%array(i), grid%lonCell%array(i), &
+ grid%sphere_radius)**2.0)
+ theta_pert = dTheta * s * sin((2.0_RKIND * pii * 0.5*(zgrid(k,i)+zgrid(k+1,i))) / L_z)
+ ! diag % theta % array(k,i) = diag % theta % array(k,i) + theta_pert
+ t(k,i) = t(k,i) + theta_pert
+ end do
+ end if
+
+
+
+ end do
+
+ !
+ ! Initialize wind field
+ !
+ allocate(psiVertex(grid % nVertices))
+ do iVtx=1,grid % nVertices
+ psiVertex(iVtx) = -grid % sphere_radius * um * ( &
+ sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
+ cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
+ )
+ end do
+ do iEdge=1,grid % nEdges
+ cell1 = grid % CellsOnEdge % array(1,iEdge)
+ cell2 = grid % CellsOnEdge % array(2,iEdge)
+ usurf = -1.0 * ( &
+ psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
+ psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
+ ) / grid%dvEdge%array(iEdge)
+ do k=1,nz1
+ ztemp = .25*( zgrid(k,cell1)+zgrid(k+1,cell1 ) &
+ +zgrid(k,cell2)+zgrid(k+1,cell2))
+
+! Top of shear layer set at 10 km
+! if(ztemp.lt.10000.) then
+ u(k,iEdge) = usurf * sqrt(1.+2.*shear*ztemp)
+! else
+! u(k,iEdge) = usurf * sqrt(1.+2.*shear*10000.)
+! end if
+ end do
+ end do
+ deallocate(psiVertex)
+
+ do k=1,nz1
+ ztemp = .5*( zw(k)+zw(k+1))
+! if(ztemp.lt.10000.) then
+ grid % u_init % array(k) = um * sqrt(1.+2.*shear*ztemp)
+! else
+! grid % u_init % array(k) = um * sqrt(1.+2.*shear*10000.)
+! end if
+ end do
+
+!
+! reference sounding based on dry atmosphere
+!
+ do i=1, grid % nCells
+
+ tsurf = t0*exp(-shear*um**2/gravity*sin(grid%latCell%array(i))**2)
+
+!! JBK fix 20120801
+!! pis = exp(-um**2*sin(grid%latCell%array(i))**2/(2.*cp*tsurf))
+ pis = 1.
+
+ pitop(i) = pis-.5*dzw(1)*gravity/(cp*tb(1,1)*zz(1,1))
+ do k=2,nz1
+ pitop(i) = pitop(i)-dzu(k)*gravity/(cp*(fzm(k)*tb(k,1)+fzp(k)*tb(k-1,1)) &
+ *(fzm(k)*zz(k,1)+fzp(k)*zz(k-1,1)))
+ end do
+ pitop(i) = pitop(i)-.5*dzw(nz1)*gravity/(cp*tb(nz1,1)*zz(nz1,1))
+ ptopb(i) = p0*pitop(i)**(1./rcp)
+
+ pb(nz1,i) = pitop(i)+.5*dzw(nz1)*gravity/(cp*tb(nz1,i)*zz(nz1,i))
+ p (nz1,i) = pitop(i)+.5*dzw(nz1)*gravity/(cp*t (nz1,i)*zz(nz1,i))
+ do k=nz1-1,1,-1
+ pb(k,i) = pb(k+1,i) + dzu(k+1)*gravity/(cp*.5*(tb(k,i)+tb(k+1,i)) &
+ *.5*(zz(k,i)+zz(k+1,i)))
+ p (k,i) = p (k+1,i) + dzu(k+1)*gravity/(cp*.5*(t (k,i)+t (k+1,i)) &
+ *.5*(zz(k,i)+zz(k+1,i)))
+ end do
+ do k=1,nz1
+ rb (k,i) = pb(k,i)**(1./rcv)/((rgas/p0)*tb(k,i)*zz(k,i))
+ rtb(k,i) = rb(k,i)*tb(k,i)
+ rr (k,i) = p (k,i)**(1./rcv)/((rgas/p0)*t (k,i)*zz(k,i))-rb(k,i)
+ cqw(k,i) = 1.
+ end do
+ end do
+
+ write(0,*) ' ***** base state sounding ***** '
+ write(0,*) 'k pb p rb rtb rr tb t'
+ do k=1,grid%nVertLevels
+ write(0,'(i2,7(2x,f14.9))') k,pb(k,1),p(k,1),rb(k,1),rtb(k,1),rr(k,1),tb(k,1),t(k,1)
+ end do
+
+ scalars(index_qv,:,:) = 0.
+!!!
+!-------------------------------------------------------------------
+! ITERATIONS TO CONVERGE MOIST SOUNDING
+ do itr=1,30
+
+ do i = 1, grid % nCells
+
+ tsurf = t0*exp(-shear*um**2/gravity*sin(grid%latCell%array(i))**2)
+ pis = exp(-um**2*sin(grid%latCell%array(i))**2/(2.*cp*tsurf))
+! pis = 1.
+
+ pitop(i) = pis-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1))
+
+ do k=2,nz1
+ pitop(i) = pitop(i)-dzu(k)*gravity/(cp*cqw(k,1)*(fzm(k)*t (k,1)+fzp(k)*t (k-1,1)) &
+ *(fzm(k)*zz(k,1)+fzp(k)*zz(k-1,1)))
+ end do
+ pitop(i) = pitop(i) - .5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,1))/(cp*t(nz1,1)*zz(nz1,1))
+ ptop = p0*pitop(i)**(1./rcp)
+
+ pp(nz1,i) = ptop-ptopb(i)+.5*dzw(nz1)*gravity* &
+ (rr(nz1,i)+(rr(nz1,i)+rb(nz1,i))*scalars(index_qv,nz1,i))
+ do k=nz1-1,1,-1
+ pp(k,i) = pp(k+1,i)+dzu(k+1)*gravity* &
+ (fzm(k)*(rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i)) &
+ +fzp(k)*(rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i)))
+ end do
+ do k=1,nz1
+ rt(k,i) = (pp(k,i)/(rgas*zz(k,i)) &
+ -rtb(k,i)*(p(k,i)-pb(k,i)))/p(k,i)
+ p (k,i) = (zz(k,i)*(rgas/p0)*(rtb(k,i)+rt(k,i)))**rcv
+ rr(k,i) = (rt(k,i)-rb(k,i)*(t(k,i)-tb(k,i)))/t(k,i)
+ end do
+!
+! update water vapor mixing ratio from humitidty profile
+!
+ do k=1,nz1
+ temp = p(k,i)*t(k,i)
+ pres = p0*p(k,i)**(1./rcp)
+ qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
+ scalars(index_qv,k,i) = min(0.014_RKIND,rh(k,i)*qvs)
+ end do
+
+ do k=1,nz1
+ t (k,i) = t(k,i)*(1.+1.61*scalars(index_qv,k,i))
+ end do
+ do k=2,nz1
+ cqw(k,i) = 1./(1.+.5*( scalars(index_qv,k-1,i) &
+ +scalars(index_qv,k ,i)))
+ end do
+
+ end do ! loop over cells
+
+ end do ! iteration loop
+!----------------------------------------------------------------------
+!
+ write(0,*) ' *** sounding for the simulation ***'
+ write(0,*) ' z theta pres qv rho_m u rr'
+ do k=1,nz1
+ write(0,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
+ t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
+ .01*p0*p(k,1)**(1./rcp), &
+ 1000.*scalars(index_qv,k,1), &
+ (rb(k,1)+rr(k,1))*(1.+scalars(index_qv,k,1)), &
+ grid % u_init % array(k), rr(k,1)
+ end do
+
+ do i=1,grid % ncells
+ do k=1,nz1
+ rho_zz(k,i) = rb(k,i)+rr(k,i)
+ end do
+
+ do k=1,nz1
+ grid % t_init % array(k,i) = t(k,i)
+ end do
+ end do
+
+ do i=1,grid % nEdges
+ cell1 = grid % CellsOnEdge % array(1,i)
+ cell2 = grid % CellsOnEdge % array(2,i)
+ if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+ do k=1,nz1
+ ru (k,i) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*u(k,i)
+ end do
+ end if
+ end do
+
+!
+! pre-calculation z-metric terms in omega eqn.
+!
+ do iEdge = 1,grid % nEdges
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+
+ do k = 1, grid%nVertLevels
+
+ if (config_theta_adv_order == 2) then
+!! test for metric consistency - forces 2nd order metrics with 4th order advection
+! if (config_theta_adv_order == 4) then
+
+ z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
+
+ else !theta_adv_order == 3 or 4
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
+ - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+ if (config_theta_adv_order == 3) then
+ z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ else
+ z_edge3 = 0.
+ end if
+
+ end if
+
+ zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
+ zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
+ zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
+ zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
+
+! if (k /= 1) then
+! zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb(k,1,iEdge)
+! zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb(k,2,iEdge)
+! zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb3(k,1,iEdge)
+! zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb3(k,2,iEdge)
+! end if
+
+ end do
+
+ end if
+ end do
+
+! for including terrain
+ state % w % array(:,:) = 0.0
+ diag % rw % array(:,:) = 0.0
+
+!
+! calculation of omega, rw = zx * ru + zz * rw
+!
+
+! do iEdge = 1,grid % nEdges
+
+! cell1 = CellsOnEdge(1,iEdge)
+! cell2 = CellsOnEdge(2,iEdge)
+
+! if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+! do k = 2, grid%nVertLevels
+! flux = (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))
+! diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + zf(k,2,iEdge)*flux
+! diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - zf(k,1,iEdge)*flux
+
+! if (config_theta_adv_order ==3) then
+! diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
+! - sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+! diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
+! + sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+! end if
+
+! end do
+! end if
+
+! end do
+
+ ! Compute w from rho_zz and rw
+ do iCell=1,grid%nCells
+ do k=2,grid%nVertLevels
+ state % w % array(k,iCell) = diag % rw % array(k,iCell) &
+ / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
+ end do
+ end do
+
+
+ do iEdge=1,grid % nEdges
+ grid % fEdge % array(iEdge) = 0.
+ end do
+
+ do iVtx=1,grid % nVertices
+ grid % fVertex % array(iVtx) = 0.
+ end do
+
+ !
+ ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
+ !
+ diag % v % array(:,:) = 0.0
+ do iEdge = 1, grid%nEdges
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ if (eoe > 0) then
+ do k = 1, grid%nVertLevels
+ diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
+ end do
+ end if
+ end do
+ end do
+
+! do k=1,grid%nVertLevels
+! write(0,*) ' k,u_init, t_init, qv_init ',k,grid % u_init % array(k),grid % t_init% array(k),grid % qv_init % array(k)
+! end do
+
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
+ end do
+ end do
+
+! MGD FOR 3-1:
+! zt = 10000.0
+! nVertLevels = 10
+! X = 125
+! dt = 12.
+! nso = 8
+! 2nd-order horiz mixing = 50.0
+
+ end subroutine init_atm_test_case_reduced_radius
+
+
+!--------------------- TEST CASE 9 -----------------------------------------------
+
+
+ subroutine init_atm_test_case_resting_atmosphere(grid, state, diag, test_case)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Setup resting atmosphere test case with terrian
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+ type (state_type), intent(inout) :: state
+ type (diag_type), intent(inout) :: diag
+ integer, intent(in) :: test_case
+
+ real (kind=RKIND), parameter :: t0=300., alpha=0.
+ real (kind=RKIND) :: hm
+
+ real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
+ real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
+ real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars, deriv_two, zb, zb3
+
+ !This is temporary variable here. It just need when calculate tangential velocity v.
+ integer :: eoe, j
+ integer, dimension(:), pointer :: nEdgesOnEdge
+ integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge, cellsOnCell, edgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcedge, AreaCell, xCell, yCell
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
+
+ integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
+ integer :: index_qv
+
+ real (kind=RKIND) :: ptop, p0, pis, flux, d2fdx2_cell1, d2fdx2_cell2
+ real(kind=RKIND), dimension(:), pointer :: hs, hs1
+
+ real (kind=RKIND) :: ztemp, zd, zt, dz, str, zh, hmax
+
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh
+ real (kind=RKIND) :: es, qvs, xnutr, ptemp
+ integer :: iter, nsm, kz
+
+ type (field1DReal), pointer :: tempField
+ type (field1DReal), target :: tempFieldTarget
+
+ type (block_type), pointer :: block
+ type (parallel_info), pointer :: parinfo
+ type (dm_info), pointer :: dminfo
+
+ real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
+ real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex
+
+ real (kind=RKIND) :: d1, d2, d3, cof1, cof2, cf1, cf2, cf3
+ real (kind=RKIND) :: um, us, rcp, rcv, gamma, xa, zinb, zint, tinv, th_inb, th_int
+ real (kind=RKIND) :: xmid, temp, pres, a_scale, rad, shear, tsurf, usurf, sm0, dzmina, dzmina_global, dzminf
+
+ real (kind=RKIND) :: xi, yi, r1m, r2m, xc, yc, xla, zinv, xn2, xn2m, xn2l, sm, dzh, dzht, dzmin, z_edge, z_edge3
+
+ integer, dimension(grid % nCells, 2) :: next_cell
+ real (kind=RKIND), dimension(grid % nCells) :: pitop, ptopb
+ logical, parameter :: hybrid = .false.
+! logical, parameter :: hybrid = .true.
+
+ block => grid % block
+ parinfo => block % parinfo
+ dminfo => block % domain % dminfo
+
+
+ !
+ ! Scale all distances
+ !
+ a_scale = grid % sphere_radius
+
+ grid % xCell % array = grid % xCell % array * a_scale
+ grid % yCell % array = grid % yCell % array * a_scale
+ grid % zCell % array = grid % zCell % array * a_scale
+ grid % xVertex % array = grid % xVertex % array * a_scale
+ grid % yVertex % array = grid % yVertex % array * a_scale
+ grid % zVertex % array = grid % zVertex % array * a_scale
+ grid % xEdge % array = grid % xEdge % array * a_scale
+ grid % yEdge % array = grid % yEdge % array * a_scale
+ grid % zEdge % array = grid % zEdge % array * a_scale
+ grid % dvEdge % array = grid % dvEdge % array * a_scale
+ grid % dcEdge % array = grid % dcEdge % array * a_scale
+ grid % areaCell % array = grid % areaCell % array * a_scale**2.0
+ grid % areaTriangle % array = grid % areaTriangle % array * a_scale**2.0
+ grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a_scale**2.0
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ AreaCell => grid % AreaCell % array
+ CellsOnEdge => grid % CellsOnEdge % array
+ cellsOnCell => grid % cellsOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ deriv_two => grid % deriv_two % array
+
+ nz1 = grid % nVertLevels
+ nz = nz1 + 1
+ nCellsSolve = grid % nCellsSolve
+
+ zgrid => grid % zgrid % array
+ zb => grid % zb % array
+ zb3 => grid % zb3 % array
+ rdzw => grid % rdzw % array
+ dzu => grid % dzu % array
+ rdzu => grid % rdzu % array
+ fzm => grid % fzm % array
+ fzp => grid % fzp % array
+ zx => grid % zx % array
+ zz => grid % zz % array
+ hx => grid % hx % array
+ dss => grid % dss % array
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+
+ ppb => diag % pressure_base % array
+ pb => diag % exner_base % array
+ rb => diag % rho_base % array
+ tb => diag % theta_base % array
+ rtb => diag % rtheta_base % array
+ p => diag % exner % array
+ cqw => diag % cqw % array
+
+ rho_zz => state % rho_zz % array
+
+ pp => diag % pressure_p % array
+ rr => diag % rho_p % array
+ t => state % theta_m % array
+ rt => diag % rtheta_p % array
+ u => state % u % array
+ ru => diag % ru % array
+
+ scalars => state % scalars % array
+
+ index_qv = state % index_qv
+
+ scalars(:,:,:) = 0.
+
+ call atm_initialize_advection_rk(grid)
+ call atm_initialize_deformation_weights(grid)
+
+ xnutr = 0.1
+ zd = 12000.
+
+ p0 = 1.e+05
+ rcp = rgas/cp
+ rcv = rgas/(cp-rgas)
+
+ ! metrics for hybrid coordinate and vertical stretching
+ str = 1.0
+
+ zt = 12000.
+
+ dz = zt/float(nz1)
+! write(0,*) ' dz = ',dz
+
+ do k=1,nz
+ zw(k) = (real(k-1)/real(nz1))**str*zt
+ if(k.gt.1) dzw(k-1) = zw(k)-zw(k-1)
+ end do
+
+! ah(k) governs the transition between terrain-following
+! and pure height coordinates
+! ah(k) = 1 is a smoothed terrain-following coordinate
+! ah(k) = 1.-zw(k)/zt is the basic terrain-following coordinate
+! ah(k) = 0 is a height coordinate
+
+ write(6,*) ' hybrid = ',hybrid
+ kz = nz
+
+ if(hybrid) then
+
+ zh = zt
+
+ do k=1,nz
+ if(zw(k).lt.zh) then
+
+! if(k.le.2) then
+! ah(k) = 1.
+! else
+! ah(k) = cos(.5*pii*(zw(k)-zw(2))/zh)**6
+! end if
+
+! ah(k) = cos(.5*pii*zw(k)/zh)**6
+ ah(k) = cos(.5*pii*zw(k)/zh)**2
+!
+! ah(k) = ah(k)*(1.-zw(k)/zt)
+!
+ else
+ ah(k) = 0.
+ kz = min(kz,k)
+ end if
+ end do
+
+ else
+
+ do k=1,nz
+ ah(k) = 1.-zw(k)/zt
+ end do
+
+ end if
+
+
+ do k=1,nz
+ write(6,*) k,zw(k), ah(k)
+ end do
+
+ write(0,*) 'EARTH RADIUS = ', grid % sphere_radius
+
+! MGD 2-0-0, not used in 2-0-1
+ if (trim(config_dcmip_case) == '2-0-0') then
+ ! for hx computation
+ r1m = .75*pii
+ r2m = pii/16.
+ end if
+
+! MGD 2-0-1, not used in 2-0-0
+ if (trim(config_dcmip_case) == '2-0-1') then
+! setting for terrain
+! xa = pii/16. ! for specifying mtn with in degrees
+ xa = pii*grid%sphere_radius/16. ! corresponds to ~11 grid intervals across entire mtn with 2 deg res
+ end if
+
+
+! MGD both 2-0-0 and 2-0-1
+ hm = 2000.0
+
+ do iCell=1,grid % nCells
+
+
+ if (trim(config_dcmip_case) == '2-0-0') then
+! Comb mountain as specified for DCMIP case 2.0
+! MGD BEGIN 2-0-0
+ xi = grid % lonCell % array(iCell)
+ yi = grid % latCell % array(iCell)
+
+ rad = acos(cos(xi)*cos(yi))
+
+ if (rad.lt.r1m) THEN
+ hx(1,iCell) = hm*cos(.5*pii*rad/r1m)**2.*cos(pii*rad/r2m)**2
+ else
+ hx(1,iCell) = 0.
+ end if
+! MGD END 2-0-0
+ end if
+
+ if (trim(config_dcmip_case) == '2-0-1') then
+! cosine**2 ridge
+! MGD BEGIN 2-0-1
+
+ xi = grid % lonCell % array(iCell)
+ yi = grid % latCell % array(iCell)
+ xc = sphere_distance(yi, xi, yi, 0., grid % sphere_radius)
+ yc = sphere_distance(yi, xi, 0., xi, grid % sphere_radius)
+
+ if (abs(xc).ge.xa) then ! for mtn ridge with uniform width in km
+! if (abs(xi).ge.xa.and.abs(2.*pii-xi).ge.xa) then ! for mtn ridge with uniform width in degrees
+ hx(1,iCell) = 0.
+ else
+! for mtn ridge with uniform width in km
+ hx(1,iCell) = hm*cos(.5*pii*xc/xa)**2*cos(yc/grid % sphere_radius)
+! for mtn ridge with uniform width in degrees
+! hx(1,iCell) = hm*cos(.5*pii*xi/xa)**2*cos(yc/grid % sphere_radius)
+ end if
+! MGD END 2-0-1
+ end if
+
+ hx(:,iCell) = hx(1,iCell)
+
+ hx(nz,iCell) = zt
+
+ end do
+
+ hmax = maxval(hx(1,:))
+ write(6,*) "max terrain height = ",hmax
+
+ if (config_smooth_surfaces) then
+
+ write(0,*) ' '
+ write(0,*) ' Smoothing vertical coordinate surfaces'
+ write(0,*) ' '
+
+ allocate(hs (grid % nCells+1))
+ allocate(hs1(grid % nCells+1))
+
+ dzmin = 0.5
+ sm0 = 0.5
+ nsm = 30
+
+ write(6,*) 'dzmin = ',dzmin,' sm0 = ',sm0,' nsm = ',nsm
+
+ do k=2,kz-1
+ hx(k,:) = hx(k-1,:)
+ dzminf = zw(k)-zw(k-1)
+
+! dzmin = max(0.5_RKIND,1.-.5*zw(k)/hm)
+
+ sm = sm0*max( min(.5*zw(k)/hm,1.0_RKIND), .05 )
+
+ do i=1,nsm
+ do iCell=1,grid % nCells
+ hs1(iCell) = 0.
+ do j = 1,nEdgesOnCell(iCell)
+
+ hs1(iCell) = hs1(iCell) + dvEdge(edgesOnCell(j,iCell)) &
+ / dcEdge(edgesOnCell(j,iCell)) &
+ * (hx(k,cellsOnCell(j,iCell))-hx(k,iCell))
+ end do
+ hs1(iCell) = hx(k,iCell) + sm*hs1(iCell)
+
+ hs(iCell) = 0.
+ ! do j = 1,nEdgesOnCell(iCell)
+ ! hs(iCell) = hs(iCell) + dvEdge(edgesOnCell(j,iCell)) &
+ ! / dcEdge(edgesOnCell(j,iCell)) &
+ ! * (hs1(cellsOnCell(j,iCell))-hs1(iCell))
+ ! end do
+ hs(iCell) = hs1(iCell) - 0.*hs(iCell)
+
+ end do
+
+ tempField => tempFieldTarget
+ tempField % block => block
+ tempField % dimSizes(1) = grid % nCells
+ tempField % sendList => parinfo % cellsToSend
+ tempField % recvList => parinfo % cellsToRecv
+ tempField % copyList => parinfo % cellsToCopy
+ tempField % array => hs
+ tempField % prev => null()
+ tempField % next => null()
+
+ call mpas_dmpar_exch_halo_field(tempField)
+
+ ! dzmina = minval(hs(:)-hx(k-1,:))
+ dzmina = minval(zw(k)+ah(k)*hs(1:grid%nCellsSolve)-zw(k-1)-ah(k-1)*hx(k-1,1:grid%nCellsSolve))
+ call mpas_dmpar_min_real(dminfo, dzmina, dzmina_global)
+ ! write(0,*) ' k,i, dzmina, dzmin, zw(k)-zw(k-1) ', k,i, dzmina, dzmin, zw(k)-zw(k-1)
+ if (dzmina_global >= dzmin*(zw(k)-zw(k-1))) then
+ hx(k,:)=hs(:)
+ dzminf = dzmina_global
+ else
+ exit
+ end if
+ end do
+ write(0,*) k,i,sm,dzminf/(zw(k)-zw(k-1)),dzmina/(zw(k)-zw(k-1))
+ end do
+
+ do k=kz,nz
+ hx(k,:) = 0.
+ end do
+
+ deallocate(hs )
+ deallocate(hs1)
+
+ else
+
+ do k=2,nz1
+ dzmina = minval(zw(k)+ah(k)*hx(k,:)-zw(k-1)-ah(k-1)*hx(k-1,:))
+ write(0,*) k,dzmina/(zw(k)-zw(k-1))
+ end do
+
+ end if
+
+
+ do iCell=1,grid % nCells
+ do k=1,nz
+ zgrid(k,iCell) = zw(k) + ah(k)*hx(k,iCell)
+ end do
+ do k=1,nz1
+ zz (k,iCell) = (zw(k+1)-zw(k))/(zgrid(k+1,iCell)-zgrid(k,iCell))
+ end do
+ end do
+
+ do i=1, grid % nEdges
+ iCell1 = grid % CellsOnEdge % array(1,i)
+ iCell2 = grid % CellsOnEdge % array(2,i)
+ do k=1,nz
+ zx (k,i) = (zgrid(k,iCell2)-zgrid(k,iCell1)) / grid % dcEdge % array(i)
+ end do
+ end do
+ do i=1, grid % nCells
+ do k=1,nz1
+ ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
+ dss(k,i) = 0.
+ ztemp = zgrid(k,i)
+ if(ztemp.gt.zd+.1) then
+ dss(k,i) = dss(k,i)+xnutr*sin(.5*pii*(ztemp-zd)/(zt-zd))**2
+ end if
+ end do
+ enddo
+
+ write(0,*) ' grid metrics setup complete '
+
+ do k=1,nz1
+ dzw (k) = zw(k+1)-zw(k)
+ rdzw(k) = 1./dzw(k)
+ zu(k ) = .5*(zw(k)+zw(k+1))
+ end do
+ do k=2,nz1
+ dzu (k) = .5*(dzw(k)+dzw(k-1))
+ rdzu(k) = 1./dzu(k)
+ fzp (k) = .5* dzw(k )/dzu(k)
+ fzm (k) = .5* dzw(k-1)/dzu(k)
+ rdzwp(k) = dzw(k-1)/(dzw(k )*(dzw(k)+dzw(k-1)))
+ rdzwm(k) = dzw(k )/(dzw(k-1)*(dzw(k)+dzw(k-1)))
+ end do
+
+! d1 = .5*dzw(1)
+! d2 = dzw(1)+.5*dzw(2)
+! d3 = dzw(1)+dzw(2)+.5*dzw(3)
+! cf1 = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+! cf2 = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+! cf3 = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+
+ cof1 = (2.*dzu(2)+dzu(3))/(dzu(2)+dzu(3))*dzw(1)/dzu(2)
+ cof2 = dzu(2) /(dzu(2)+dzu(3))*dzw(1)/dzu(3)
+ cf1 = fzp(2) + cof1
+ cf2 = fzm(2) - cof1 - cof2
+ cf3 = cof2
+
+ grid % cf1 % scalar = cf1
+ grid % cf2 % scalar = cf2
+ grid % cf3 % scalar = cf3
+
+ um = 0.
+ gamma = .0065 ! temp lapse rate in K/km
+
+! MGD BEGIN 2-0-0
+ if (trim(config_dcmip_case) == '2-0-0') then
+ zinb = zt ! no inversion layer
+ zint = zt ! no inversion layer
+ end if
+! MGD END 2-0-0
+! MGD BEGIN 2-0-1
+ if (trim(config_dcmip_case) == '2-0-1') then
+ zinb = 3000. ! bottom of inversion layer
+ zint = 5000. ! top of inversion layer
+ end if
+! MGD END 2-0-1
+
+ ! computing intermediate T and Theta used to build sounding that includes inversion layer
+ tinv = t0-gamma*zinb
+ th_inb = t0*(1.-gamma*zinb/t0)**(1.-gravity/(cp*gamma))
+ th_int = th_inb*exp((gravity*(zint-zinb))/(cp*tinv))
+ write(6,*) ' zinb = ',zinb,' zint = ',zint,' tinv = ',tinv,'th_inb = ',th_inb,' th_int = ',th_int
+
+ do i=1,grid % nCells
+
+ pis = 1.
+
+ do k=1,nz1
+ ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
+
+! Isothermal reference sounding
+
+ tb(k,i) = t0*exp(gravity*ztemp/(cp*t0))
+
+! Low level inversion initial sounding
+
+ if(ztemp.le.zinb) then
+ t (k,i) = t0*(1.-gamma*ztemp/t0)**(1.-gravity/(cp*gamma))
+ else if(ztemp.le.zint) then
+ t (k,i) = th_inb*exp((gravity*(ztemp-zinb))/(cp*tinv))
+ else
+ t (k,i) = th_int*(1.-gamma*(ztemp-zint)/tinv)**(1.-gravity/(cp*gamma))
+ end if
+
+ rh(k,i) = 0.
+ end do
+ end do
+
+ !
+ ! Initialize wind field
+ !
+ do iEdge=1,grid % nEdges
+ do k=1,nz1
+ u(k,iEdge) = um
+ end do
+ end do
+
+ do k=1,nz1
+ grid % u_init % array(k) = um
+ end do
+
+!
+! reference sounding based on dry atmosphere
+!
+ do i=1, grid % nCells
+
+ pis = 1.
+
+ pitop(i) = pis-.5*dzw(1)*gravity/(cp*tb(1,1)*zz(1,1))
+ do k=2,nz1
+ pitop(i) = pitop(i)-dzu(k)*gravity/(cp*(fzm(k)*tb(k,1)+fzp(k)*tb(k-1,1)) &
+ *(fzm(k)*zz(k,1)+fzp(k)*zz(k-1,1)))
+ end do
+ pitop(i) = pitop(i)-.5*dzw(nz1)*gravity/(cp*tb(nz1,1)*zz(nz1,1))
+ ptopb(i) = p0*pitop(i)**(1./rcp)
+
+ pb(nz1,i) = pitop(i)+.5*dzw(nz1)*gravity/(cp*tb(nz1,i)*zz(nz1,i))
+ p (nz1,i) = pitop(i)+.5*dzw(nz1)*gravity/(cp*t (nz1,i)*zz(nz1,i))
+ do k=nz1-1,1,-1
+ pb(k,i) = pb(k+1,i) + dzu(k+1)*gravity/(cp*.5*(tb(k,i)+tb(k+1,i)) &
+ *.5*(zz(k,i)+zz(k+1,i)))
+ p (k,i) = p (k+1,i) + dzu(k+1)*gravity/(cp*.5*(t (k,i)+t (k+1,i)) &
+ *.5*(zz(k,i)+zz(k+1,i)))
+ end do
+ do k=1,nz1
+ rb (k,i) = pb(k,i)**(1./rcv)/((rgas/p0)*tb(k,i)*zz(k,i))
+ rtb(k,i) = rb(k,i)*tb(k,i)
+ rr (k,i) = p (k,i)**(1./rcv)/((rgas/p0)*t (k,i)*zz(k,i))-rb(k,i)
+ cqw(k,i) = 1.
+ end do
+ end do
+
+ write(0,*) ' ***** base state sounding ***** '
+ write(0,*) 'k pb p rb rtb rr tb t'
+ do k=1,grid%nVertLevels
+ write(0,'(i2,7(2x,f14.9))') k,pb(k,1),p(k,1),rb(k,1),rtb(k,1),rr(k,1),tb(k,1),t(k,1)
+ end do
+
+ scalars(index_qv,:,:) = 0.
+!!!
+!-------------------------------------------------------------------
+! ITERATIONS TO CONVERGE MOIST SOUNDING
+ do itr=1,30
+
+ do i = 1, grid % nCells
+
+ pis = 1.
+
+ pitop(i) = pis-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1))
+
+ do k=2,nz1
+ pitop(i) = pitop(i)-dzu(k)*gravity/(cp*cqw(k,1)*(fzm(k)*t (k,1)+fzp(k)*t (k-1,1)) &
+ *(fzm(k)*zz(k,1)+fzp(k)*zz(k-1,1)))
+ end do
+ pitop(i) = pitop(i) - .5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,1))/(cp*t(nz1,1)*zz(nz1,1))
+ ptop = p0*pitop(i)**(1./rcp)
+
+ pp(nz1,i) = ptop-ptopb(i)+.5*dzw(nz1)*gravity* &
+ (rr(nz1,i)+(rr(nz1,i)+rb(nz1,i))*scalars(index_qv,nz1,i))
+ do k=nz1-1,1,-1
+ pp(k,i) = pp(k+1,i)+dzu(k+1)*gravity* &
+ (fzm(k)*(rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i)) &
+ +fzp(k)*(rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i)))
+ end do
+ do k=1,nz1
+ rt(k,i) = (pp(k,i)/(rgas*zz(k,i)) &
+ -rtb(k,i)*(p(k,i)-pb(k,i)))/p(k,i)
+ p (k,i) = (zz(k,i)*(rgas/p0)*(rtb(k,i)+rt(k,i)))**rcv
+ rr(k,i) = (rt(k,i)-rb(k,i)*(t(k,i)-tb(k,i)))/t(k,i)
+ end do
+!
+! update water vapor mixing ratio from humitidty profile
+!
+ do k=1,nz1
+ temp = p(k,i)*t(k,i)
+ pres = p0*p(k,i)**(1./rcp)
+ qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
+ scalars(index_qv,k,i) = min(0.014_RKIND,rh(k,i)*qvs)
+ end do
+
+ do k=1,nz1
+ t (k,i) = t(k,i)*(1.+1.61*scalars(index_qv,k,i))
+ end do
+ do k=2,nz1
+ cqw(k,i) = 1./(1.+.5*( scalars(index_qv,k-1,i) &
+ +scalars(index_qv,k ,i)))
+ end do
+
+ end do ! loop over cells
+
+ end do ! iteration loop
+!----------------------------------------------------------------------
+!
+ write(0,*) ' *** sounding for the simulation ***'
+ write(0,*) ' z temp theta pres rho_m u rr'
+ do k=1,nz1
+ write(0,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
+ t(k,1)/(1.+1.61*scalars(index_qv,k,1))*p(k,1), &
+ t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
+ .01*p0*p(k,1)**(1./rcp), &
+! 1000.*scalars(index_qv,k,1), &
+ (rb(k,1)+rr(k,1))*(1.+scalars(index_qv,k,1)), &
+ grid % u_init % array(k), rr(k,1)
+ end do
+
+ do i=1,grid % ncells
+ do k=1,nz1
+ rho_zz(k,i) = rb(k,i)+rr(k,i)
+ end do
+
+ do k=1,nz1
+ grid % t_init % array(k,i) = t(k,i)
+ end do
+ end do
+
+ do i=1,grid % nEdges
+ cell1 = grid % CellsOnEdge % array(1,i)
+ cell2 = grid % CellsOnEdge % array(2,i)
+ if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+ do k=1,nz1
+ ru (k,i) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*u(k,i)
+ end do
+ end if
+ end do
+
+!
+! pre-calculation z-metric terms in omega eqn.
+!
+ do iEdge = 1,grid % nEdges
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+
+ do k = 1, grid%nVertLevels
+
+ if (config_theta_adv_order == 2) then
+!! test for metric consistency - forces 2nd order metrics with 4th order advection
+! if (config_theta_adv_order == 4) then
+
+ z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
+
+ else !theta_adv_order == 3 or 4
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
+ - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+ if (config_theta_adv_order == 3) then
+ z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ else
+ z_edge3 = 0.
+ end if
+
+ end if
+
+ zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
+ zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
+ zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
+ zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
+
+! if (k /= 1) then
+! zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb(k,1,iEdge)
+! zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb(k,2,iEdge)
+! zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb3(k,1,iEdge)
+! zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb3(k,2,iEdge)
+! end if
+
+ end do
+
+ end if
+ end do
+
+! for including terrain
+ state % w % array(:,:) = 0.0
+ diag % rw % array(:,:) = 0.0
+
+!
+! calculation of omega, rw = zx * ru + zz * rw
+!
+
+! do iEdge = 1,grid % nEdges
+
+! cell1 = CellsOnEdge(1,iEdge)
+! cell2 = CellsOnEdge(2,iEdge)
+
+! if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+! do k = 2, grid%nVertLevels
+! flux = (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))
+! diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + zf(k,2,iEdge)*flux
+! diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - zf(k,1,iEdge)*flux
+
+! if (config_theta_adv_order ==3) then
+! diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
+! - sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+! diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
+! + sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+! end if
+
+! end do
+! end if
+
+! end do
+
+ ! Compute w from rho_zz and rw
+ do iCell=1,grid%nCells
+ do k=2,grid%nVertLevels
+ state % w % array(k,iCell) = diag % rw % array(k,iCell) &
+ / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
+ end do
+ end do
+
+
+ do iEdge=1,grid % nEdges
+ grid % fEdge % array(iEdge) = 0.
+ end do
+
+ do iVtx=1,grid % nVertices
+ grid % fVertex % array(iVtx) = 0.
+ end do
+
+ !
+ ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
+ !
+ diag % v % array(:,:) = 0.0
+ do iEdge = 1, grid%nEdges
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ if (eoe > 0) then
+ do k = 1, grid%nVertLevels
+ diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
+ end do
+ end if
+ end do
+ end do
+
+! do k=1,grid%nVertLevels
+! write(0,*) ' k,u_init, t_init, qv_init ',k,grid % u_init % array(k),grid % t_init% array(k),grid % qv_init % array(k)
+! end do
+
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
+ end do
+ end do
+
+ end subroutine init_atm_test_case_resting_atmosphere
+
+
integer function nearest_cell(target_lat, target_lon, &
start_cell, &
nCells, maxEdges, nEdgesOnCell, cellsOnCell, latCell, lonCell)
Modified: branches/omp_blocks/openmp_test/src/core_nhyd_atmos/Makefile
===================================================================
--- branches/omp_blocks/openmp_test/src/core_nhyd_atmos/Makefile 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_nhyd_atmos/Makefile 2012-09-07 15:07:43 UTC (rev 2145)
@@ -4,7 +4,6 @@
#PHYSICS=
OBJS = mpas_atm_mpas_core.o \
- mpas_atm_test_cases.o \
mpas_atm_time_integration.o \
mpas_atm_advection.o
@@ -19,13 +18,11 @@
core_hyd: $(OBJS)
ar -ru libdycore.a $(OBJS) phys/*.o
-mpas_atm_test_cases.o: mpas_atm_advection.o
-
mpas_atm_time_integration.o:
mpas_atm_advection.o:
-mpas_atm_mpas_core.o: mpas_atm_advection.o mpas_atm_test_cases.o mpas_atm_time_integration.o
+mpas_atm_mpas_core.o: mpas_atm_advection.o mpas_atm_time_integration.o
clean:
( cd ../core_atmos_physics; make clean )
Modified: branches/omp_blocks/openmp_test/src/core_nhyd_atmos/Registry
===================================================================
--- branches/omp_blocks/openmp_test/src/core_nhyd_atmos/Registry 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_nhyd_atmos/Registry 2012-09-07 15:07:43 UTC (rev 2145)
@@ -1,7 +1,6 @@
%
% namelist type namelist_record name default_value
%
-namelist integer nhyd_model config_test_case 0
namelist character nhyd_model config_time_integration SRK3
namelist real nhyd_model config_dt 600.0
namelist character nhyd_model config_calendar_type gregorian
@@ -38,7 +37,6 @@
namelist integer nhyd_model config_num_halos 2
namelist real damping config_zd 22000.0
namelist real damping config_xnutr 0.0
-namelist integer dimensions config_nvertlevels 26
namelist character io config_input_name init.nc
namelist character io config_sfc_update_name sfc_update.nc
namelist character io config_output_name output.nc
@@ -69,7 +67,7 @@
dim FIFTEEN 15
dim TWENTYONE 21
dim R3 3
-dim nVertLevels namelist:config_nvertlevels
+dim nVertLevels nVertLevels
dim nVertLevelsP1 nVertLevels+1
%
Modified: branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_advection.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_advection.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_advection.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -111,9 +111,9 @@
do i=1,n
advCells(i+1,iCell) = cell_list(i)
- xc(i) = grid % xCell % array(advCells(i+1,iCell))/a
- yc(i) = grid % yCell % array(advCells(i+1,iCell))/a
- zc(i) = grid % zCell % array(advCells(i+1,iCell))/a
+ xc(i) = grid % xCell % array(advCells(i+1,iCell))/grid%sphere_radius
+ yc(i) = grid % yCell % array(advCells(i+1,iCell))/grid%sphere_radius
+ zc(i) = grid % zCell % array(advCells(i+1,iCell))/grid%sphere_radius
end do
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
@@ -131,8 +131,8 @@
xc(i+1), yc(i+1), zc(i+1), &
xc(ip2), yc(ip2), zc(ip2) )
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1) )
+ dl_sphere(i) = grid%sphere_radius*arc_length( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1) )
end do
length_scale = 1.
@@ -262,12 +262,12 @@
if (ip1 > n-1) ip1 = 1
iEdge = grid % EdgesOnCell % array (i,iCell)
- xv1 = grid % xVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- yv1 = grid % yVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- zv1 = grid % zVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- xv2 = grid % xVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/a
+ xv1 = grid % xVertex % array(grid % verticesOnEdge % array (1,iedge))/grid%sphere_radius
+ yv1 = grid % yVertex % array(grid % verticesOnEdge % array (1,iedge))/grid%sphere_radius
+ zv1 = grid % zVertex % array(grid % verticesOnEdge % array (1,iedge))/grid%sphere_radius
+ xv2 = grid % xVertex % array(grid % verticesOnEdge % array (2,iedge))/grid%sphere_radius
+ yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/grid%sphere_radius
+ zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/grid%sphere_radius
if ( grid % on_a_sphere ) then
call arc_bisect( xv1, yv1, zv1, &
@@ -825,16 +825,16 @@
! compute poynomial fit for this cell if all needed neighbors exist
if (grid % on_a_sphere) then
- xc(1) = grid % xCell % array(iCell)/a
- yc(1) = grid % yCell % array(iCell)/a
- zc(1) = grid % zCell % array(iCell)/a
+ xc(1) = grid % xCell % array(iCell)/grid%sphere_radius
+ yc(1) = grid % yCell % array(iCell)/grid%sphere_radius
+ zc(1) = grid % zCell % array(iCell)/grid%sphere_radius
do i=2,n
iv = grid % verticesOnCell % array(i-1,iCell)
- xc(i) = grid % xVertex % array(iv)/a
- yc(i) = grid % yVertex % array(iv)/a
- zc(i) = grid % zVertex % array(iv)/a
+ xc(i) = grid % xVertex % array(iv)/grid%sphere_radius
+ yc(i) = grid % yVertex % array(iv)/grid%sphere_radius
+ zc(i) = grid % zVertex % array(iv)/grid%sphere_radius
end do
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
@@ -852,8 +852,8 @@
xc(i+1), yc(i+1), zc(i+1), &
xc(ip2), yc(ip2), zc(ip2) )
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1) )
+ dl_sphere(i) = grid%sphere_radius*arc_length( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1) )
end do
length_scale = 1.
Modified: branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_mpas_core.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_mpas_core.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_mpas_core.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -20,7 +20,6 @@
use mpas_configure
use mpas_kind_types
use mpas_grid_types
- use atm_test_cases
implicit none
@@ -37,8 +36,21 @@
integer :: i
integer :: ierr
- if (.not. config_do_restart) call atm_setup_test_case(domain)
+ if (.not. config_do_restart) then
+ ! Code that was previously handled by atm_setup_test_case()
+
+ block => domain % blocklist
+ do while (associated(block))
+ do i=2,nTimeLevs
+ call mpas_copy_state(block % state % time_levs(i) % state, block % state % time_levs(1) % state)
+ end do
+ block => block % next
+ end do
+
+ end if
+
+
!
! Initialize core
!
Deleted: branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_test_cases.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_test_cases.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_test_cases.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -1,2184 +0,0 @@
-module atm_test_cases
-
- use mpas_grid_types
- use mpas_configure
- use mpas_constants
- use mpas_dmpar
- use atm_advection
-#ifdef DO_PHYSICS
- use mpas_atmphys_control
-#endif
-
-
- contains
-
-
- subroutine atm_setup_test_case(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Configure grid metadata and model state for the hydrostatic test case
- ! specified in the namelist
- !
- ! Output: block - a subset (not necessarily proper) of the model domain to be
- ! initialized
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: i
- type (block_type), pointer :: block_ptr
-
- if (config_test_case == 0) then
- write(0,*) ' Using initial conditions from input file'
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
- block_ptr => block_ptr % next
- end do
-
- 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))
- write(0,*) ' calling test case setup '
- call atm_test_case_jw(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, &
- block_ptr % diag_physics ,config_test_case)
- write(0,*) ' returned from test case setup '
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- else if ((config_test_case == 4) .or. (config_test_case ==5)) then
-
- write(0,*) ' squall line - super cell test case '
- if (config_test_case == 4) write(0,*) ' squall line test case'
- if (config_test_case == 5) write(0,*) ' supercell test case'
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- write(0,*) ' calling test case setup '
- call atm_test_case_squall_line(domain % dminfo, block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
- write(0,*) ' returned from test case setup '
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 6 ) then
-
- write(0,*) ' mountain wave test case '
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- write(0,*) ' calling test case setup '
- call atm_test_case_mtn_wave(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
- write(0,*) ' returned from test case setup '
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- else
-
-
- write(0,*) ' Only test case 1, 2, 3, 4, 5 and 6 are currently supported for nonhydrostatic core '
- stop
- end if
-
-
-#ifdef DO_PHYSICS
- !initialization of surface input variables technically not needed to run our current set of
- !idealized test cases:
- if (config_test_case > 0) then
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call physics_idealized_init(block_ptr % mesh, block_ptr % sfc_input)
- block_ptr => block_ptr % next
- end do
-
- endif
-#endif
-
- end subroutine atm_setup_test_case
-
-!----------------------------------------------------------------------------------------------------------
-
- subroutine atm_test_case_jw(grid, state, diag, diag_physics, test_case)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
- type (diag_type), intent(inout) :: diag
- type (diag_physics_type), intent(inout) :: diag_physics
- integer, intent(in) :: test_case
-
- real (kind=RKIND), parameter :: u0 = 35.0
- real (kind=RKIND), parameter :: alpha_grid = 0. ! no grid rotation
- real (kind=RKIND), parameter :: omega_e = 7.29212e-05
- real (kind=RKIND), parameter :: t0b = 250., t0 = 288., delta_t = 4.8e+05, dtdz = 0.005, eta_t = 0.2
- 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 :: k_x = 9. ! Normal mode wave number
-
- real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
- real (kind=RKIND), dimension(:), pointer :: surface_pressure
- real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx
- real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt
- real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
-
-!.. initialization of moisture:
- integer:: index_qv
- real (kind=RKIND),parameter :: rh_max = 0.40 ! Maximum relative humidity
-! real (kind=RKIND),parameter :: rh_max = 0.70 ! Maximum relative humidity
- real (kind=RKIND),dimension(:,:), pointer:: qsat, relhum
- real (kind=RKIND),dimension(:,:,:),pointer:: scalars
-!.. end initialization of moisture.
-
- integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
-
- !This is temporary variable here. It just need when calculate tangential velocity v.
- integer :: eoe, j
- integer, dimension(:), pointer :: nEdgesOnEdge
- integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge
- real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
-
- 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) :: r_earth, etavs, ztemp, zd, zt, dz, gam, delt, str
-
- real (kind=RKIND) :: es, qvs, xnutr, znut, ptemp
- integer :: iter
-
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: znuc, znuv, bn, divh, dpn
-
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: sh, zw, ah
- real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
- real (kind=RKIND), dimension(grid % nVertLevels ) :: eta, etav, teta, ppi, tt, temperature_1d
-
- real (kind=RKIND) :: d1, d2, d3, cf1, cf2, cf3, cof1, cof2, psurf
-
- ! storage for (lat,z) arrays for zonal velocity calculation
-
- logical, parameter :: rebalance = .true.
- integer, parameter :: nlat=721
- real (kind=RKIND), dimension(grid % nVertLevels) :: flux_zonal
- real (kind=RKIND), dimension(grid % nVertLevels + 1, nlat) :: zgrid_2d
- real (kind=RKIND), dimension(grid % nVertLevels, nlat) :: u_2d, pp_2d, rho_2d, qv_2d, etavs_2d, zz_2d
- real (kind=RKIND), dimension(grid % nVertLevels, nlat-1) :: zx_2d
- real (kind=RKIND), dimension(nlat) :: lat_2d
- real (kind=RKIND) :: dlat, hx_1d
- real (kind=RKIND) :: z_edge, z_edge3, d2fdx2_cell1, d2fdx2_cell2
-
- logical, parameter :: moisture = .true.
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
- weightsOnEdge => grid % weightsOnEdge % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- dvEdge => grid % dvEdge % array
- AreaCell => grid % AreaCell % array
- CellsOnEdge => grid % CellsOnEdge % array
-
- deriv_two => grid % deriv_two % array
- zb => grid % zb % array
- zb3 => grid % zb3% array
-
- nz1 = grid % nVertLevels
- nz = nz1 + 1
- nCellsSolve = grid % nCellsSolve
-
- zgrid => grid % zgrid % array
- rdzw => grid % rdzw % array
- dzu => grid % dzu % array
- rdzu => grid % rdzu % array
- fzm => grid % fzm % array
- fzp => grid % fzp % array
- zx => grid % zx % array
- zz => grid % zz % array
- hx => grid % hx % array
- dss => grid % dss % array
-
- pb => diag % exner_base % array
- rb => diag % rho_base % array
- tb => diag % theta_base % array
- rtb => diag % rtheta_base % array
- p => diag % exner % array
-
- ppb => diag % pressure_base % array
- pp => diag % pressure_p % array
-
- rho_zz => state % rho_zz % array
- rr => diag % rho_p % array
- t => state % theta_m % array
- rt => diag % rtheta_p % array
-
- surface_pressure => diag % surface_pressure % array
-
-!.. initialization of moisture:
- scalars => state % scalars % array
- qsat => diag_physics % qsat % array
- relhum => diag_physics % relhum % array
- scalars(:,:,:) = 0.0
- qsat(:,:) = 0.0
- relhum(:,:) = 0.0
- qv_2d(:,:) = 0.0
-!.. end initialization of moisture.
-
- surface_pressure(:) = 0.0
-
- call atm_initialize_advection_rk(grid)
- call atm_initialize_deformation_weights(grid)
-
- index_qv = state % index_qv
-
- xnutr = 0.
- zd = 12000.
- znut = eta_t
-
- etavs = (1.-0.252)*pii/2.
- r_earth = a
- p0 = 1.e+05
-
- write(0,*) ' point 1 in test case setup '
-
-! We may pass in an hx(:,:) that has been precomputed elsewhere.
-! For now it is independent of k
-
- do iCell=1,grid % nCells
- do k=1,nz
- phi = grid % latCell % array (iCell)
- hx(k,iCell) = u0/gravity*cos(etavs)**1.5 &
- *((-2.*sin(phi)**6 &
- *(cos(phi)**2+1./3.)+10./63.) &
- *(u0)*cos(etavs)**1.5 &
- +(1.6*cos(phi)**3 &
- *(sin(phi)**2+2./3.)-pii/4.)*r_earth*omega_e)
- enddo
- enddo
-
- ! Metrics for hybrid coordinate and vertical stretching
-
- str = 1.8
- zt = 45000.
- dz = zt/float(nz1)
-
- write(0,*) ' hx computation complete '
-
- do k=1,nz
-
-! sh(k) is the stretching specified for height surfaces
-
- sh(k) = (real(k-1)*dz/zt)**str
-
-! to specify specific heights zc(k) for coordinate surfaces,
-! input zc(k) and define sh(k) = zc(k)/zt
-! zw(k) is the hieght of zeta surfaces
-! zw(k) = (k-1)*dz yields constant dzeta
-! and nonconstant dzeta/dz
-! zw(k) = sh(k)*zt yields nonconstant dzeta
-! and nearly constant dzeta/dz
-
- zw(k) = float(k-1)*dz
-! zw(k) = sh(k)*zt
-!
-! ah(k) governs the transition between terrain-following
-! and pureheight coordinates
-! ah(k) = 0 is a terrain-following coordinate
-! ah(k) = 1 is a height coordinate
-
- ah(k) = 1.-cos(.5*pii*(k-1)*dz/zt)**6
-! ah(k) = 0.
- write(0,*) ' k, sh, zw, ah ',k,sh(k),zw(k),ah(k)
- end do
- do k=1,nz1
- dzw (k) = zw(k+1)-zw(k)
- rdzw(k) = 1./dzw(k)
- zu(k ) = .5*(zw(k)+zw(k+1))
- end do
- do k=2,nz1
- dzu (k) = .5*(dzw(k)+dzw(k-1))
- rdzu(k) = 1./dzu(k)
- fzp (k) = .5* dzw(k )/dzu(k)
- fzm (k) = .5* dzw(k-1)/dzu(k)
- rdzwp(k) = dzw(k-1)/(dzw(k )*(dzw(k)+dzw(k-1)))
- rdzwm(k) = dzw(k )/(dzw(k-1)*(dzw(k)+dzw(k-1)))
- end do
-
-!********** how are we storing cf1, cf2 and cf3?
-
- COF1 = (2.*DZU(2)+DZU(3))/(DZU(2)+DZU(3))*DZW(1)/DZU(2)
- COF2 = DZU(2) /(DZU(2)+DZU(3))*DZW(1)/DZU(3)
- CF1 = FZP(2) + COF1
- CF2 = FZM(2) - COF1 - COF2
- CF3 = COF2
-
-! d1 = .5*dzw(1)
-! d2 = dzw(1)+.5*dzw(2)
-! d3 = dzw(1)+dzw(2)+.5*dzw(3)
-! cf1 = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
-! cf2 = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
-! cf3 = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
-
- write(0,*) ' cf1, cf2, cf3 = ',cf1,cf2,cf3
-
- grid % cf1 % scalar = cf1
- grid % cf2 % scalar = cf2
- grid % cf3 % scalar = cf3
-
- do iCell=1,grid % nCells
- do k=1,nz
- zgrid(k,iCell) = (1.-ah(k))*(sh(k)*(zt-hx(k,iCell))+hx(k,iCell)) &
- + ah(k) * sh(k)* zt
- end do
- do k=1,nz1
- zz (k,iCell) = (zw(k+1)-zw(k))/(zgrid(k+1,iCell)-zgrid(k,iCell))
- end do
- end do
-
- do i=1, grid % nEdges
- iCell1 = grid % CellsOnEdge % array(1,i)
- iCell2 = grid % CellsOnEdge % array(2,i)
- do k=1,nz
- zx (k,i) = (zgrid(k,iCell2)-zgrid(k,iCell1)) / grid % dcEdge % array(i)
- end do
- end do
- do i=1, grid % nCells
- do k=1,nz1
- ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
- dss(k,i) = 0.
- ztemp = zgrid(k,i)
- if(ztemp.gt.zd+.1) then
- dss(k,i) = dss(k,i)+xnutr*sin(.5*pii*(ztemp-zd)/(zt-zd))**2
- end if
- end do
- enddo
-
- !do k=1,nz1
- ! write(0,*) ' k, zgrid(k,1),hx(k,1) ',k,zgrid(k,1),hx(k,1)
- !enddo
-
- !do k=1,nz1
- ! write(0,*) ' k, zx(k,1) ',k,zx(k,1)
- !enddo
-
- write(0,*) ' grid metrics setup complete '
-
-!************** section for 2d (z,lat) calc for zonal velocity
-
- dlat = 0.5*pii/float(nlat-1)
- do i = 1,nlat
-
- lat_2d(i) = float(i-1)*dlat
- phi = lat_2d(i)
- hx_1d = u0/gravity*cos(etavs)**1.5 &
- *((-2.*sin(phi)**6 &
- *(cos(phi)**2+1./3.)+10./63.) &
- *(u0)*cos(etavs)**1.5 &
- +(1.6*cos(phi)**3 &
- *(sin(phi)**2+2./3.)-pii/4.)*r_earth*omega_e)
-
- do k=1,nz
- zgrid_2d(k,i) = (1.-ah(k))*(sh(k)*(zt-hx_1d)+hx_1d) &
- + ah(k) * sh(k)* zt
- end do
- do k=1,nz1
- zz_2d (k,i) = (zw(k+1)-zw(k))/(zgrid_2d(k+1,i)-zgrid_2d(k,i))
- end do
-
- do k=1,nz1
- ztemp = .5*(zgrid_2d(k+1,i)+zgrid_2d(k,i))
- ppb(k,i) = p0*exp(-gravity*ztemp/(rgas*t0b))
- pb (k,i) = (ppb(k,i)/p0)**(rgas/cp)
- rb (k,i) = ppb(k,i)/(rgas*t0b*zz_2d(k,i))
- tb (k,i) = t0b/pb(k,i)
- rtb(k,i) = rb(k,i)*tb(k,i)
- p (k,i) = pb(k,i)
- pp (k,i) = 0.
- rr (k,i) = 0.
- end do
-
-
- do itr = 1,10
-
- do k=1,nz1
- eta (k) = (ppb(k,i)+pp(k,i))/p0
- etav(k) = (eta(k)-.252)*pii/2.
- if(eta(k).ge.znut) then
- teta(k) = t0*eta(k)**(rgas*dtdz/gravity)
- else
- teta(k) = t0*eta(k)**(rgas*dtdz/gravity) + delta_t*(znut-eta(k))**5
- end if
- end do
-
- phi = lat_2d (i)
- do k=1,nz1
- temperature_1d(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
- *sqrt(cos(etav(k)))* &
- ((-2.*sin(phi)**6 &
- *(cos(phi)**2+1./3.)+10./63.) &
- *2.*u0*cos(etav(k))**1.5 &
- +(1.6*cos(phi)**3 &
- *(sin(phi)**2+2./3.)-pii/4.)*r_earth*omega_e)/(1.+0.61*qv_2d(k,i))
-
-
- ztemp = .5*(zgrid_2d(k,i)+zgrid_2d(k+1,i))
- ptemp = ppb(k,i) + pp(k,i)
-
- !get moisture
- if (moisture) then
- qv_2d(k,i) = env_qv( ztemp, temperature_1d(k), ptemp, rh_max )
- end if
-
- tt(k) = temperature_1d(k)*(1.+1.61*qv_2d(k,i))
- end do
-
- do itrp = 1,25
- do k=1,nz1
- rr(k,i) = (pp(k,i)/(rgas*zz_2d(k,i)) - rb(k,i)*(tt(k)-t0b))/tt(k)
- end do
-
- ppi(1) = p0-.5*dzw(1)*gravity &
- *(1.25*(rr(1,i)+rb(1,i))*(1.+qv_2d(1,i)) &
- -.25*(rr(2,i)+rb(2,i))*(1.+qv_2d(2,i)))
-
- ppi(1) = ppi(1)-ppb(1,i)
- do k=1,nz1-1
- ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity* &
- (rr(k ,i)+(rr(k ,i)+rb(k ,i))*qv_2d(k ,i) &
- +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*qv_2d(k+1,i))
- end do
-
- do k=1,nz1
- pp(k,i) = .2*ppi(k)+.8*pp(k,i)
- end do
-
- end do ! end inner iteration loop itrp
-
- end do ! end outer iteration loop itr
-
- do k=1,nz1
- rho_2d(k,i) = rr(k,i)+rb(k,i)
- pp_2d(k,i) = pp(k,i)
- etavs_2d(k,i) = ((ppb(k,i)+pp(k,i))/p0 - 0.252)*pii/2.
- u_2d(k,i) = u0*(sin(2.*lat_2d(i))**2) *(cos(etavs_2d(k,i))**1.5)
- end do
-
- end do ! end loop over latitudes for 2D zonal wind field calc
-
- !SHP-balance:: in case of rebalacing for geostrophic wind component
- if (rebalance) then
-
- do i=1,nlat-1
- do k=1,nz1
- zx_2d(k,i) = (zgrid_2d(k,i+1)-zgrid_2d(k,i))/(dlat*r_earth)
- end do
- end do
-
- call atm_recompute_geostrophic_wind(u_2d,rho_2d,pp_2d,qv_2d,lat_2d,zz_2d,zx_2d, &
- cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat)
-
- end if
-
-!******************************************************************
-
-!
-!---- baroclinc wave initialization ---------------------------------
-!
-! reference sounding based on dry isothermal atmosphere
-!
- do i=1, grid % nCells
- do k=1,nz1
- ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
- ppb(k,i) = p0*exp(-gravity*ztemp/(rgas*t0b))
- pb (k,i) = (ppb(k,i)/p0)**(rgas/cp)
- rb (k,i) = ppb(k,i)/(rgas*t0b*zz(k,i))
- tb (k,i) = t0b/pb(k,i)
- rtb(k,i) = rb(k,i)*tb(k,i)
- p (k,i) = pb(k,i)
- pp (k,i) = 0.
- rr (k,i) = 0.
- end do
-
-! if(i == 1) then
-! do k=1,nz1
-! write(0,*) ' k, ppb, pb, rb, tb (k,1) ',k,ppb(k,1),pb(k,1),rb(k,1)*zz(k,1),tb(k,1)
-! enddo
-! end if
-
- 200 format(4i6,8(1x,e15.8))
- 201 format(3i6,8(1x,e15.8))
- 202 format(2i6,10(1x,e15.8))
- 203 format(i6,10(1x,e15.8))
-
-! iterations to converge temperature as a function of pressure
-!
- do itr = 1,10
-
- do k=1,nz1
- eta (k) = (ppb(k,i)+pp(k,i))/p0
- etav(k) = (eta(k)-.252)*pii/2.
- if(eta(k).ge.znut) then
- teta(k) = t0*eta(k)**(rgas*dtdz/gravity)
- else
- teta(k) = t0*eta(k)**(rgas*dtdz/gravity) + delta_t*(znut-eta(k))**5
- end if
- end do
- phi = grid % latCell % array (i)
- do k=1,nz1
- temperature_1d(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
- *sqrt(cos(etav(k)))* &
- ((-2.*sin(phi)**6 &
- *(cos(phi)**2+1./3.)+10./63.) &
- *2.*u0*cos(etav(k))**1.5 &
- +(1.6*cos(phi)**3 &
- *(sin(phi)**2+2./3.)-pii/4.)*r_earth*omega_e)/(1.+0.61*scalars(index_qv,k,i))
-
- ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
- ptemp = ppb(k,i) + pp(k,i)
-
- !get moisture
- if (moisture) then
-
- !scalars(index_qv,k,i) = env_qv( ztemp, temperature_1d(k), ptemp, rh_max )
-
- if(ptemp < 50000.) then
- relhum(k,i) = 0.0
- elseif(ptemp > p0) then
- relhum(k,i) = 1.0
- else
- relhum(k,i) = (1.-((p0-ptemp)/50000.)**1.25)
- endif
- relhum(k,i) = min(rh_max,relhum(k,i))
-
- !.. calculation of water vapor mixing ratio:
- if (temperature_1d(k) > 273.15) then
- es = 1000.*0.6112*exp(17.67*(temperature_1d(k)-273.15)/(temperature_1d(k)-29.65))
- else
- es = 1000.*0.6112*exp(21.8745584*(temperature_1d(k)-273.15)/(temperature_1d(k)-7.66))
- end if
- qsat(k,i) = (287.04/461.6)*es/(ptemp-es)
- if(relhum(k,i) .eq. 0.0) qsat(k,i) = 0.0
- scalars(index_qv,k,i) = relhum(k,i)*qsat(k,i)
- end if
-
- tt(k) = temperature_1d(k)*(1.+1.61*scalars(index_qv,k,i))
-
- end do
-
- do itrp = 1,25
- do k=1,nz1
- rr(k,i) = (pp(k,i)/(rgas*zz(k,i)) - rb(k,i)*(tt(k)-t0b))/tt(k)
- end do
-
- ppi(1) = p0-.5*dzw(1)*gravity &
- *(1.25*(rr(1,i)+rb(1,i))*(1.+scalars(index_qv,1,i)) &
- -.25*(rr(2,i)+rb(2,i))*(1.+scalars(index_qv,2,i)))
-
- ppi(1) = ppi(1)-ppb(1,i)
- do k=1,nz1-1
- ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity* &
- (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i) &
- +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))
- end do
-
- do k=1,nz1
- pp(k,i) = .2*ppi(k)+.8*pp(k,i)
- end do
-
- end do ! end inner iteration loop itrp
-
- end do ! end outer iteration loop itr
-
- do k=1,nz1
- p (k,i) = ((ppb(k,i)+pp(k,i))/p0)**(rgas/cp)
- t (k,i) = tt(k)/p(k,i)
- rt (k,i) = t(k,i)*rr(k,i)+rb(k,i)*(t(k,i)-tb(k,i))
- rho_zz (k,i) = rb(k,i) + rr(k,i)
- end do
-
- !calculation of surface pressure:
- surface_pressure(i) = 0.5*dzw(1)*gravity &
- * (1.25*(rr(1,i) + rb(1,i)) * (1. + scalars(index_qv,1,i)) &
- - 0.25*(rr(2,i) + rb(2,i)) * (1. + scalars(index_qv,2,i)))
- surface_pressure(i) = surface_pressure(i) + pp(1,i) + ppb(1,i)
-
- end do ! end loop over cells
-
- !write(0,*)
- !write(0,*) '--- initialization of water vapor:'
- !do iCell = 1, grid % nCells
- ! if(iCell == 1 .or. iCell == grid % nCells) then
- ! do k = nz1, 1, -1
- ! write(0,202) iCell,k,t(k,iCell),relhum(k,iCell),qsat(k,iCell),scalars(index_qv,k,iCell)
- ! enddo
- ! write(0,*)
- ! endif
- !enddo
-
- lat_pert = latitude_pert*pii/180.
- lon_pert = longitude_pert*pii/180.
-
- do iEdge=1,grid % nEdges
-
- vtx1 = grid % VerticesOnEdge % array (1,iEdge)
- vtx2 = grid % VerticesOnEdge % array (2,iEdge)
- lat1 = grid%latVertex%array(vtx1)
- lat2 = grid%latVertex%array(vtx2)
- iCell1 = grid % cellsOnEdge % array(1,iEdge)
- iCell2 = grid % cellsOnEdge % array(2,iEdge)
- flux = (0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1)))*a/grid % dvEdge % array(iEdge)
-
- if (config_test_case == 2) then
- r_pert = sphere_distance( grid % latEdge % array (iEdge), grid % lonEdge % array (iEdge), &
- lat_pert, lon_pert, 1.0_RKIND)/(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
-
- if (rebalance) then
-
- call atm_calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1,lat2,grid % dvEdge % array(iEdge),a,u0,nz1,nlat)
- do k=1,grid % nVertLevels
- fluxk = u0*flux_zonal(k)/(0.5*(rb(k,iCell1)+rb(k,iCell2)+rr(k,iCell1)+rr(k,iCell2)))
- state % u % array(k,iEdge) = fluxk + u_pert
- end do
-
- else
-
- do k=1,grid % nVertLevels
- etavs = (0.5*(ppb(k,iCell1)+ppb(k,iCell2)+pp(k,iCell1)+pp(k,iCell2))/p0 - 0.252)*pii/2.
- fluxk = u0*flux*(cos(etavs)**1.5)
- state % u % array(k,iEdge) = fluxk + u_pert
- end do
-
- end if
-
- cell1 = grid % CellsOnEdge % array(1,iEdge)
- cell2 = grid % CellsOnEdge % array(2,iEdge)
- do k=1,nz1
- diag % ru % array (k,iEdge) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*state % u % array (k,iEdge)
- end do
-
- !
- ! Generate rotated Coriolis field
- !
-
- grid % fEdge % array(iEdge) = 2.0 * omega * &
- ( -cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha_grid) + &
- sin(grid%latEdge%array(iEdge)) * cos(alpha_grid) &
- )
- end do
-
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 2.0 * omega * &
- (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha_grid) + &
- sin(grid%latVertex%array(iVtx)) * cos(alpha_grid) &
- )
- end do
-
- !
- ! CALCULATION OF OMEGA, RW = ZX * RU + ZZ * RW
- !
-
- !
- ! pre-calculation z-metric terms in omega eqn.
- !
- do iEdge = 1,grid % nEdges
- cell1 = CellsOnEdge(1,iEdge)
- cell2 = CellsOnEdge(2,iEdge)
-
- do k = 1, grid%nVertLevels
-
- if (config_theta_adv_order == 2) then
-
- z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
-
- else if (config_theta_adv_order == 3 .or. config_theta_adv_order ==4) then !theta_adv_order == 3 or 4
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
- d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
- - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- if (config_theta_adv_order == 3) then
- z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
- else
- z_edge3 = 0.
- end if
-
- end if
-
- zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
- zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
- zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
- zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
-
- end do
-
- end do
-
- ! for including terrain
- diag % rw % array = 0.
- state % w % array = 0.
- do iEdge = 1,grid % nEdges
-
- cell1 = CellsOnEdge(1,iEdge)
- cell2 = CellsOnEdge(2,iEdge)
-
- do k = 2, grid%nVertLevels
- flux = (fzm(k)*diag % ru % array(k,iEdge)+fzp(k)*diag % ru % array(k-1,iEdge))
- diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + (fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb(k,2,iEdge)*flux
- diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - (fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb(k,1,iEdge)*flux
-
- if (config_theta_adv_order ==3) then
- diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
- - sign(1.0_RKIND,diag % ru % array(k,iEdge))*config_coef_3rd_order* &
- (fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb3(k,2,iEdge)*flux
- diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
- + sign(1.0_RKIND,diag % ru % array(k,iEdge))*config_coef_3rd_order* &
- (fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb3(k,1,iEdge)*flux
- end if
-
- end do
-
- end do
-
- ! Compute w from rho_zz and rw
- do iCell=1,grid%nCells
- do k=2,grid%nVertLevels
- state % w % array(k,iCell) = diag % rw % array(k,iCell) &
- / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
- end do
- end do
-
-
- !
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
- diag % v % array(:,:) = 0.0
- do iEdge = 1, grid%nEdges
- do i=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(i,iEdge)
- do k = 1, grid%nVertLevels
- diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
- end do
- end do
- end do
-
- do i=1,10
- psurf = (cf1*(ppb(1,i)+pp(1,i)) + cf2*(ppb(2,i)+pp(2,i)) + cf3*(ppb(3,i)+pp(3,i)))/100.
-
- psurf = (ppb(1,i)+pp(1,i)) + .5*dzw(1)*gravity &
- *(1.25*(rr(1,i)+rb(1,i))*(1.+scalars(index_qv,1,i)) &
- -.25*(rr(2,i)+rb(2,i))*(1.+scalars(index_qv,2,i)))
-
- write(0,*) ' i, psurf, lat ',i,psurf,grid%latCell%array(i)*180./3.1415828
- enddo
-
- ! Compute rho and theta from rho_zz and theta_m
- do iCell=1,grid%nCells
- do k=1,grid%nVertLevels
- diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
- diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
- end do
- end do
-
- end subroutine atm_test_case_jw
-
- subroutine atm_calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1_in,lat2_in,dvEdge,a,u0,nz1,nlat)
-
- implicit none
- integer, intent(in) :: nz1,nlat
- real (kind=RKIND), dimension(nz1,nlat), intent(in) :: u_2d,etavs_2d
- real (kind=RKIND), dimension(nlat), intent(in) :: lat_2d
- real (kind=RKIND), dimension(nz1), intent(out) :: flux_zonal
- real (kind=RKIND), intent(in) :: lat1_in, lat2_in, dvEdge, a, u0
-
- integer :: k,i
- real (kind=RKIND) :: lat1, lat2, w1, w2
- real (kind=RKIND) :: dlat,da,db
-
- lat1 = abs(lat1_in)
- lat2 = abs(lat2_in)
- if(lat2 <= lat1) then
- lat1 = abs(lat2_in)
- lat2 = abs(lat1_in)
- end if
-
- do k=1,nz1
- flux_zonal(k) = 0.
- end do
-
- do i=1,nlat-1
- if( (lat1 <= lat_2d(i+1)) .and. (lat2 >= lat_2d(i)) ) then
-
- dlat = lat_2d(i+1)-lat_2d(i)
- da = (max(lat1,lat_2d(i))-lat_2d(i))/dlat
- db = (min(lat2,lat_2d(i+1))-lat_2d(i))/dlat
- w1 = (db-da) -0.5*(db-da)**2
- w2 = 0.5*(db-da)**2
-
- do k=1,nz1
- flux_zonal(k) = flux_zonal(k) + w1*u_2d(k,i) + w2*u_2d(k,i+1)
- end do
-
- end if
-
- end do
-
-! renormalize for setting cell-face fluxes
-
- do k=1,nz1
- flux_zonal(k) = sign(1.0_RKIND,lat2_in-lat1_in)*flux_zonal(k)*dlat*a/dvEdge/u0
- end do
-
- end subroutine atm_calc_flux_zonal
-
- !SHP-balance
- subroutine atm_recompute_geostrophic_wind(u_2d,rho_2d,pp_2d,qv_2d,lat_2d,zz_2d,zx_2d, &
- cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat)
-
- implicit none
- integer, intent(in) :: nz1,nlat
- real (kind=RKIND), dimension(nz1,nlat), intent(inout) :: u_2d
- real (kind=RKIND), dimension(nz1,nlat), intent(in) :: rho_2d, pp_2d, qv_2d, zz_2d
- real (kind=RKIND), dimension(nz1,nlat-1), intent(in) :: zx_2d
- real (kind=RKIND), dimension(nlat), intent(in) :: lat_2d
- real (kind=RKIND), dimension(nz1), intent(in) :: fzm, fzp, rdzw
- real (kind=RKIND), intent(in) :: cf1, cf2, cf3, dlat
-
- !local variable
- real (kind=RKIND), dimension(nz1,nlat-1) :: pgrad, ru, u
- real (kind=RKIND), dimension(nlat-1) :: f
- real (kind=RKIND), dimension(nz1+1) :: dpzx
-
- real (kind=RKIND), parameter :: omega_e = 7.29212e-05
- real (kind=RKIND) :: rdx, qtot, r_earth, phi
- integer :: k,i, itr
-
- r_earth = a
- rdx = 1./(dlat*r_earth)
-
- do i=1,nlat-1
- do k=1,nz1
- pgrad(k,i) = rdx*(pp_2d(k,i+1)/zz_2d(k,i+1)-pp_2d(k,i)/zz_2d(k,i))
- end do
-
- dpzx(:) = 0.
-
- k=1
- dpzx(k) = .5*zx_2d(k,i)*(cf1*(pp_2d(k ,i+1)+pp_2d(k ,i)) &
- +cf2*(pp_2d(k+1,i+1)+pp_2d(k+1,i)) &
- +cf3*(pp_2d(k+2,i+1)+pp_2d(k+2,i)))
- do k=2,nz1
- dpzx(k) = .5*zx_2d(k,i)*(fzm(k)*(pp_2d(k ,i+1)+pp_2d(k ,i)) &
- +fzp(k)*(pp_2d(k-1,i+1)+pp_2d(k-1,i)))
- end do
-
- do k=1,nz1
- pgrad(k,i) = pgrad(k,i) - rdzw(k)*(dpzx(k+1)-dpzx(k))
- end do
- end do
-
-
- !initial value of v and rv -> that is from analytic sln.
- do i=1,nlat-1
- do k=1,nz1
- u(k,i) = .5*(u_2d(k,i)+u_2d(k,i+1))
- ru(k,i) = u(k,i)*(rho_2d(k,i)+rho_2d(k,i+1))*.5
- end do
- end do
-
- write(0,*) "MAX U wind before REBALANCING ---->", maxval(abs(u))
-
- !re-calculate geostrophic wind using iteration
- do itr=1,50
- do i=1,nlat-1
- phi = (lat_2d(i)+lat_2d(i+1))/2.
- f(i) = 2.*omega_e*sin(phi)
- do k=1,nz1
- if (f(i).eq.0.) then
- ru(k,i) = 0.
- else
- qtot = .5*(qv_2d(k,i)+qv_2d(k,i+1))
- ru(k,i) = - ( 1./(1.+qtot)*pgrad(k,i) + tan(phi)/r_earth*u(k,i)*ru(k,i) )/f(i)
- end if
- u(k,i) = ru(k,i)*2./(rho_2d(k,i)+rho_2d(k,i+1))
- end do
- end do
- end do
-
- write(0,*) "MAX U wind after REBALANCING ---->", maxval(abs(u))
-
- !update 2d ru
- do i=2,nlat-1
- do k=1,nz1
- u_2d(k,i) = (ru(k,i-1)+ru(k,i))*.5
- end do
- end do
-
- i=1
- do k=1,nz1
- u_2d(k,i) = (3.*u_2d(k,i+1)-u_2d(k,i+2))*.5
- end do
- i=nlat
- do k=1,nz1
- u_2d(k,i) = (3.*u_2d(k,i-1)-u_2d(k,i-2))*.5
- end do
-
-
- end subroutine atm_recompute_geostrophic_wind
-!----------------------------------------------------------------------------------------------------------
-
- subroutine atm_test_case_squall_line(dminfo, grid, state, diag, test_case)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup squall line and supercell test case
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (dm_info), intent(in) :: dminfo
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
- type (diag_type), intent(inout) :: diag
- integer, intent(in) :: test_case
-
- real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
- real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
- real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
- real (kind=RKIND), dimension(:,:,:), pointer :: scalars
-
- !This is temporary variable here. It just need when calculate tangential velocity v.
- integer :: eoe, j
- integer, dimension(:), pointer :: nEdgesOnEdge
- integer, dimension(:,:), pointer :: edgesOnEdge
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
-
- integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, cell1, cell2, nCellsSolve
- integer :: index_qv
-
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: znu, znw, znwc, znwv
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: znuc, znuv
-
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
- real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
-
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh, thi, tbi, cqwb
-
- real (kind=RKIND) :: r, xnutr
- real (kind=RKIND) :: ztemp, zd, zt, dz, str
-
- real (kind=RKIND), dimension(grid % nVertLevels ) :: qvb
- real (kind=RKIND), dimension(grid % nVertLevels ) :: t_init_1d
-
- real (kind=RKIND) :: d1, d2, d3, cf1, cf2, cf3, cof1, cof2
- real (kind=RKIND) :: ztr, thetar, ttr, thetas, um, us, zts, pitop, pibtop, ptopb, ptop, rcp, rcv, p0
- real (kind=RKIND) :: radx, radz, zcent, xmid, delt, xloc, rad, yloc, ymid, a_scale
- real (kind=RKIND) :: pres, temp, es, qvs
-
- !
- ! Scale all distances
- !
-
- a_scale = 1.0
-
- grid % xCell % array = grid % xCell % array * a_scale
- grid % yCell % array = grid % yCell % array * a_scale
- grid % zCell % array = grid % zCell % array * a_scale
- grid % xVertex % array = grid % xVertex % array * a_scale
- grid % yVertex % array = grid % yVertex % array * a_scale
- grid % zVertex % array = grid % zVertex % array * a_scale
- grid % xEdge % array = grid % xEdge % array * a_scale
- grid % yEdge % array = grid % yEdge % array * a_scale
- grid % zEdge % array = grid % zEdge % array * a_scale
- grid % dvEdge % array = grid % dvEdge % array * a_scale
- grid % dcEdge % array = grid % dcEdge % array * a_scale
- grid % areaCell % array = grid % areaCell % array * a_scale**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a_scale**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a_scale**2.0
-
- weightsOnEdge => grid % weightsOnEdge % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
-
- nz1 = grid % nVertLevels
- nz = nz1 + 1
- nCellsSolve = grid % nCellsSolve
-
- zgrid => grid % zgrid % array
- rdzw => grid % rdzw % array
- dzu => grid % dzu % array
- rdzu => grid % rdzu % array
- fzm => grid % fzm % array
- fzp => grid % fzp % array
- zx => grid % zx % array
- zz => grid % zz % array
- hx => grid % hx % array
- dss => grid % dss % array
-
- ppb => diag % pressure_base % array
- pb => diag % exner_base % array
- rb => diag % rho_base % array
- tb => diag % theta_base % array
- rtb => diag % rtheta_base % array
- p => diag % exner % array
- cqw => diag % cqw % array
-
- rho_zz => state % rho_zz % array
-
- pp => diag % pressure_p % array
- rr => diag % rho_p % array
- t => state % theta_m % array
- rt => diag % rtheta_p % array
- u => state % u % array
- ru => diag % ru % array
-
- scalars => state % scalars % array
-
- index_qv = state % index_qv
-
- scalars(:,:,:) = 0.
-
- call atm_initialize_advection_rk(grid)
- call atm_initialize_deformation_weights(grid)
-
- xnutr = 0.
- zd = 12000.
-
- p0 = 1.e+05
- rcp = rgas/cp
- rcv = rgas/(cp-rgas)
-
- write(0,*) ' point 1 in test case setup '
-
-! We may pass in an hx(:,:) that has been precomputed elsewhere.
-! For now it is independent of k
-
- do iCell=1,grid % nCells
- do k=1,nz
- hx(k,iCell) = 0. ! squall line or supercell on flat plane
- enddo
- enddo
-
- ! metrics for hybrid coordinate and vertical stretching
-
- str = 1.0
- zt = 20000.
- dz = zt/float(nz1)
-
-! write(0,*) ' dz = ',dz
- write(0,*) ' hx computation complete '
-
- do k=1,nz
-
-! sh(k) is the stretching specified for height surfaces
-
- zc(k) = zt*(real(k-1)*dz/zt)**str
-
-! to specify specific heights zc(k) for coordinate surfaces,
-! input zc(k)
-! zw(k) is the hieght of zeta surfaces
-! zw(k) = (k-1)*dz yields constant dzeta
-! and nonconstant dzeta/dz
-! zw(k) = sh(k)*zt yields nonconstant dzeta
-! and nearly constant dzeta/dz
-
-! zw(k) = float(k-1)*dz
- zw(k) = zc(k)
-!
-! ah(k) governs the transition between terrain-following
-! and pureheight coordinates
-! ah(k) = 0 is a terrain-following coordinate
-! ah(k) = 1 is a height coordinate
-
-! ah(k) = 1.-cos(.5*pii*(k-1)*dz/zt)**6
- ah(k) = 1.
-! write(0,*) ' k, zc, zw, ah ',k,zc(k),zw(k),ah(k)
- end do
- do k=1,nz1
- dzw (k) = zw(k+1)-zw(k)
- rdzw(k) = 1./dzw(k)
- zu(k ) = .5*(zw(k)+zw(k+1))
- end do
- do k=2,nz1
- dzu (k) = .5*(dzw(k)+dzw(k-1))
- rdzu(k) = 1./dzu(k)
- fzp (k) = .5* dzw(k )/dzu(k)
- fzm (k) = .5* dzw(k-1)/dzu(k)
- rdzwp(k) = dzw(k-1)/(dzw(k )*(dzw(k)+dzw(k-1)))
- rdzwm(k) = dzw(k )/(dzw(k-1)*(dzw(k)+dzw(k-1)))
- end do
-
-!********** how are we storing cf1, cf2 and cf3?
-
- COF1 = (2.*DZU(2)+DZU(3))/(DZU(2)+DZU(3))*DZW(1)/DZU(2)
- COF2 = DZU(2) /(DZU(2)+DZU(3))*DZW(1)/DZU(3)
- CF1 = FZP(2) + COF1
- CF2 = FZM(2) - COF1 - COF2
- CF3 = COF2
-
-! d1 = .5*dzw(1)
-! d2 = dzw(1)+.5*dzw(2)
-! d3 = dzw(1)+dzw(2)+.5*dzw(3)
-! cf1 = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
-! cf2 = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
-! cf3 = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
-
- grid % cf1 % scalar = cf1
- grid % cf2 % scalar = cf2
- grid % cf3 % scalar = cf3
-
- do iCell=1,grid % nCells
- do k=1,nz
- zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(k,iCell)/zt)+hx(k,iCell)) &
- + (1.-ah(k)) * zc(k)
- end do
- do k=1,nz1
- zz (k,iCell) = (zw(k+1)-zw(k))/(zgrid(k+1,iCell)-zgrid(k,iCell))
- end do
- end do
-
- do i=1, grid % nEdges
- iCell1 = grid % CellsOnEdge % array(1,i)
- iCell2 = grid % CellsOnEdge % array(2,i)
- do k=1,nz
- zx (k,i) = (zgrid(k,iCell2)-zgrid(k,iCell1)) / grid % dcEdge % array(i)
- end do
- end do
- do i=1, grid % nCells
- do k=1,nz1
- ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
- dss(k,i) = 0.
- ztemp = zgrid(k,i)
- if(ztemp.gt.zd+.1) then
- dss(k,i) = dss(k,i)+xnutr*sin(.5*pii*(ztemp-zd)/(zt-zd))**2
- end if
- end do
- enddo
-
-!
-! convective initialization
-!
- ztr = 12000.
- thetar = 343.
- ttr = 213.
- thetas = 300.5
-
-! write(0,*) ' rgas, cp, gravity ',rgas,cp, gravity
-
- if ( config_test_case == 4) then ! squall line parameters
- um = 12.
- us = 10.
- zts = 2500.
- else if (config_test_case == 5) then !supercell parameters
- um = 30.
- us = 15.
- zts = 5000.
- end if
-
- do i=1,grid % nCells
- do k=1,nz1
- ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
- if(ztemp .gt. ztr) then
- t (k,i) = thetar*exp(9.8*(ztemp-ztr)/(1003.*ttr))
- rh(k,i) = 0.25
- else
- t (k,i) = 300.+43.*(ztemp/ztr)**1.25
- rh(k,i) = (1.-0.75*(ztemp/ztr)**1.25)
- if(t(k,i).lt.thetas) t(k,i) = thetas
- end if
- tb(k,i) = t(k,i)
- thi(k,i) = t(k,i)
- tbi(k,i) = t(k,i)
- cqw(k,i) = 1.
- cqwb(k,i) = 1.
- end do
- end do
-
-! rh(:,:) = 0.
-
-! set the velocity field - we are on a plane here.
-
- do i=1, grid % nEdges
- cell1 = grid % CellsOnEdge % array(1,i)
- cell2 = grid % CellsOnEdge % array(2,i)
- if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,nz1
- ztemp = .25*( zgrid(k,cell1 )+zgrid(k+1,cell1 ) &
- +zgrid(k,cell2)+zgrid(k+1,cell2))
- if(ztemp.lt.zts) then
- u(k,i) = um*ztemp/zts
- else
- u(k,i) = um
- end if
- if(i == 1 ) grid % u_init % array(k) = u(k,i) - us
- u(k,i) = cos(grid % angleEdge % array(i)) * (u(k,i) - us)
- end do
- end if
- end do
-
- call mpas_dmpar_bcast_reals(dminfo, nz1, grid % u_init % array)
-
-!
-! for reference sounding
-!
- do itr=1,30
-
- pitop = 1.-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1))
- pibtop = 1.-.5*dzw(1)*gravity*(1.+qvb(1))/(cp*tb(1,1)*zz(1,1))
- do k=2,nz1
- pitop = pitop-dzu(k)*gravity/(cp*cqw(k,1)*.5*(t(k,1)+t(k-1,1)) &
- *.5*(zz(k,1)+zz(k-1,1)))
- pibtop = pibtop-dzu(k)*gravity/(cp*cqwb(k,1)*.5*(tb(k,1)+tb(k-1,1)) &
- *.5*(zz(k,1)+zz(k-1,1)))
-
- !write(0,*) k,pitop,tb(k,1),dzu(k),tb(k,1)
- end do
- pitop = pitop-.5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,1))/(cp*t(nz1,1)*zz(nz1,1))
- pibtop = pibtop-.5*dzw(nz1)*gravity*(1.+qvb(nz1))/(cp*tb(nz1,1)*zz(nz1,1))
-
- call mpas_dmpar_bcast_real(dminfo, pitop)
- call mpas_dmpar_bcast_real(dminfo, pibtop)
-
- ptopb = p0*pibtop**(1./rcp)
- write(6,*) 'ptopb = ',.01*ptopb
-
- do i=1, grid % nCells
- pb(nz1,i) = pibtop+.5*dzw(nz1)*gravity*(1.+qvb(nz1))/(cp*tb(nz1,i)*zz(nz1,i))
- p (nz1,i) = pitop+.5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,i))/(cp*t (nz1,i)*zz(nz1,i))
- do k=nz1-1,1,-1
- pb(k,i) = pb(k+1,i) + dzu(k+1)*gravity/(cp*cqwb(k+1,i)*.5*(tb(k,i)+tb(k+1,i)) &
- *.5*(zz(k,i)+zz(k+1,i)))
- p (k,i) = p (k+1,i) + dzu(k+1)*gravity/(cp*cqw(k+1,i)*.5*(t (k,i)+t (k+1,i)) &
- *.5*(zz(k,i)+zz(k+1,i)))
- end do
- do k=1,nz1
- rb (k,i) = pb(k,i)**(1./rcv)/((rgas/p0)*tb(k,i)*zz(k,i))
- rtb(k,i) = rb(k,i)*tb(k,i)
- rr (k,i) = p (k,i)**(1./rcv)/((rgas/p0)*t (k,i)*zz(k,i))-rb(k,i)
- ppb(k,i) = p0*(zz(k,i)*rgas*rtb(k,i)/p0)**(cp/cv)
- end do
- end do
-
- !
- ! update water vapor mixing ratio from humidity profile
- !
- do i= 1,grid%nCells
- do k=1,nz1
- temp = p(k,i)*thi(k,i)
- pres = p0*p(k,i)**(1./rcp)
- qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
- scalars(index_qv,k,i) = min(0.014_RKIND,rh(k,i)*qvs)
- end do
- end do
-
- do k=1,nz1
-!*********************************************************************
-! QVB = QV INCLUDES MOISTURE IN REFERENCE STATE
-! qvb(k) = scalars(index_qv,k,1)
-
-! QVB = 0 PRODUCES DRY REFERENCE STATE
- qvb(k) = 0.
-!*********************************************************************
- end do
-
- do i= 1,grid%nCells
- do k=1,nz1
- t (k,i) = thi(k,i)*(1.+1.61*scalars(index_qv,k,i))
- tb(k,i) = tbi(k,i)*(1.+1.61*qvb(k))
- end do
- do k=2,nz1
- cqw (k,i) = 1./(1.+.5*(scalars(index_qv,k,i)+scalars(index_qv,k-1,i)))
- cqwb(k,i) = 1./(1.+.5*(qvb(k)+qvb(k-1)))
- end do
- end do
-
- end do !end of iteration loop
-
- write(0,*) ' base state sounding '
- write(0,*) ' k, pb, rb, tb, rtb, t, rr, p, qvb'
- do k=1,grid%nVertLevels
- write (0,'(i2,8(2x,f19.15))') k,pb(k,1),rb(k,1),tb(k,1),rtb(k,1),t(k,1),rr(k,1),p(k,1),qvb(k)
- end do
-
-!
-! potential temperature perturbation
-!
-! delt = -10.
-! delt = -0.01
- delt = 3.
- radx = 10000.
- radz = 1500.
- zcent = 1500.
-
- if (config_test_case == 4) then ! squall line prameters
- call mpas_dmpar_max_real(dminfo, maxval(grid % xCell % array(:)), xmid)
- xmid = xmid * 0.5
- ymid = 0.0 ! Not used for squall line
- else if (config_test_case == 5) then ! supercell parameters
- call mpas_dmpar_max_real(dminfo, maxval(grid % xCell % array(:)), xmid)
- call mpas_dmpar_max_real(dminfo, maxval(grid % yCell % array(:)), ymid)
- xmid = xmid * 0.5
- ymid = ymid * 0.5
- end if
-
- do i=1, grid % nCells
- xloc = grid % xCell % array(i) - xmid
- if (config_test_case == 4) then
- yloc = 0. !squall line setting
- else if (config_test_case == 5) then
- yloc = grid % yCell % array(i) - ymid !supercell setting
- end if
-
- do k = 1,nz1
- ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
- rad =sqrt((xloc/radx)**2+(yloc/radx)**2+((ztemp-zcent)/radz)**2)
- if(rad.lt.1) then
- thi(k,i) = thi(k,i) + delt*cos(.5*pii*rad)**2
- end if
- t (k,i) = thi(k,i)*(1.+1.61*scalars(index_qv,k,i))
- end do
- end do
-
- do itr=1,30
-
- pitop = 1.-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1))
- do k=2,nz1
- pitop = pitop-dzu(k)*gravity/(cp*cqw(k,1)*.5*(t (k,1)+t (k-1,1)) &
- *.5*(zz(k,1)+zz(k-1,1)))
- end do
- pitop = pitop - .5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,1))/(cp*t(nz1,1)*zz(nz1,1))
- ptop = p0*pitop**(1./rcp)
- write(0,*) 'ptop = ',.01*ptop, .01*ptopb
-
- call mpas_dmpar_bcast_real(dminfo, ptop)
-
- do i = 1, grid % nCells
-
- pp(nz1,i) = ptop-ptopb+.5*dzw(nz1)*gravity* &
- (rr(nz1,i)+(rr(nz1,i)+rb(nz1,i))*scalars(index_qv,nz1,i))
- do k=nz1-1,1,-1
-! pp(k,i) = pp(k+1,i)+.5*dzu(k+1)*gravity* &
-! (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i) &
-! +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))
- pp(k,i) = pp(k+1,i)+dzu(k+1)*gravity*( &
- fzm(k+1)*(rb(k+1,i)*(scalars(index_qv,k+1,i)-qvb(k+1)) &
- +rr(k+1,i)*(1.+scalars(index_qv,k+1,i))) &
- +fzp(k+1)*(rb(k ,i)*(scalars(index_qv,k ,i)-qvb(k)) &
- +rr(k ,i)*(1.+scalars(index_qv,k ,i))))
- end do
- if (itr==1.and.i==1) then
- do k=1,nz1
- write(0,*) "pp-check", pp(k,i)
- end do
- end if
- do k=1,nz1
- rt(k,i) = (pp(k,i)/(rgas*zz(k,i)) &
- -rtb(k,i)*(p(k,i)-pb(k,i)))/p(k,i)
- p (k,i) = (zz(k,i)*(rgas/p0)*(rtb(k,i)+rt(k,i)))**rcv
- rr(k,i) = (rt(k,i)-rb(k,i)*(t(k,i)-tb(k,i)))/t(k,i)
- end do
-
- end do ! loop over cells
-
- end do ! iteration loop
-!----------------------------------------------------------------------
-!
- do k=1,nz1
- grid % qv_init % array(k) = scalars(index_qv,k,1)
- end do
-
- t_init_1d(:) = t(:,1)
- call mpas_dmpar_bcast_reals(dminfo, nz1, t_init_1d)
- call mpas_dmpar_bcast_reals(dminfo, nz1, grid % qv_init % array)
-
- do i=1,grid % ncells
- do k=1,nz1
- grid % t_init % array(k,i) = t_init_1d(k)
- rho_zz(k,i) = rb(k,i)+rr(k,i)
- end do
- end do
-
- do i=1,grid % nEdges
- cell1 = grid % CellsOnEdge % array(1,i)
- cell2 = grid % CellsOnEdge % array(2,i)
- if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,nz1
- ru (k,i) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*u(k,i)
- end do
- end if
- end do
-
-
- !
- ! we are assuming w and rw are zero for this initialization
- ! i.e., no terrain
- !
- diag % rw % array = 0.
- state % w % array = 0.
-
- grid % zb % array = 0.
- grid % zb3% array = 0.
-
- !
- ! Generate rotated Coriolis field
- !
- do iEdge=1,grid % nEdges
- grid % fEdge % array(iEdge) = 0.
- end do
-
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 0.
- end do
-
- !
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
- diag % v % array(:,:) = 0.0
- do iEdge = 1, grid%nEdges
- do i=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
- do k = 1, grid%nVertLevels
- diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
- end do
- end if
- end do
- end do
-
- ! write(0,*) ' k,u_init, t_init, qv_init '
- ! do k=1,grid%nVertLevels
- ! write(0,'(i2,3(2x,f14.10)') k,grid % u_init % array(k),grid % t_init% array(k),grid % qv_init % array(k)
- ! end do
-
- ! Compute rho and theta from rho_zz and theta_m
- do iCell=1,grid%nCells
- do k=1,grid%nVertLevels
- diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
- diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
- end do
- end do
-
- end subroutine atm_test_case_squall_line
-
-
-!----------------------------------------------------------------------------------------------------------
-
-
- subroutine atm_test_case_mtn_wave(grid, state, diag, test_case)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
- type (diag_type), intent(inout) :: diag
- integer, intent(in) :: test_case
-
- real (kind=RKIND), parameter :: t0=288., hm=250.
-
- real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
- real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
- real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
- real (kind=RKIND), dimension(:,:,:), pointer :: scalars, deriv_two, zb, zb3
-
- !This is temporary variable here. It just need when calculate tangential velocity v.
- integer :: eoe, j
- integer, dimension(:), pointer :: nEdgesOnEdge
- integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge
- real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
-
- integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
- integer :: index_qv
-
- real (kind=RKIND) :: ptop, pitop, ptopb, p0, flux, d2fdx2_cell1, d2fdx2_cell2
-
- real (kind=RKIND) :: ztemp, zd, zt, dz, str
-
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh
- real (kind=RKIND) :: es, qvs, xnutr, ptemp
- integer :: iter
-
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
- real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
-
- real (kind=RKIND) :: d1, d2, d3, cof1, cof2, cf1, cf2, cf3
- real (kind=RKIND) :: um, us, rcp, rcv
- real (kind=RKIND) :: xmid, temp, pres, a_scale
-
- real (kind=RKIND) :: xi, xa, xc, xla, zinv, xn2, xn2m, xn2l, sm, dzh, dzht, dzmin, z_edge, z_edge3
-
- integer, dimension(grid % nCells, 2) :: next_cell
- real (kind=RKIND), dimension(grid % nCells) :: hxzt
- logical, parameter :: terrain_smooth = .false.
-
- !
- ! Scale all distances
- !
-
- a_scale = 1.0
-
- grid % xCell % array = grid % xCell % array * a_scale
- grid % yCell % array = grid % yCell % array * a_scale
- grid % zCell % array = grid % zCell % array * a_scale
- grid % xVertex % array = grid % xVertex % array * a_scale
- grid % yVertex % array = grid % yVertex % array * a_scale
- grid % zVertex % array = grid % zVertex % array * a_scale
- grid % xEdge % array = grid % xEdge % array * a_scale
- grid % yEdge % array = grid % yEdge % array * a_scale
- grid % zEdge % array = grid % zEdge % array * a_scale
- grid % dvEdge % array = grid % dvEdge % array * a_scale
- grid % dcEdge % array = grid % dcEdge % array * a_scale
- grid % areaCell % array = grid % areaCell % array * a_scale**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a_scale**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a_scale**2.0
-
- weightsOnEdge => grid % weightsOnEdge % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- dvEdge => grid % dvEdge % array
- AreaCell => grid % AreaCell % array
- CellsOnEdge => grid % CellsOnEdge % array
- deriv_two => grid % deriv_two % array
-
- nz1 = grid % nVertLevels
- nz = nz1 + 1
- nCellsSolve = grid % nCellsSolve
-
- zgrid => grid % zgrid % array
- zb => grid % zb % array
- zb3 => grid % zb3 % array
- rdzw => grid % rdzw % array
- dzu => grid % dzu % array
- rdzu => grid % rdzu % array
- fzm => grid % fzm % array
- fzp => grid % fzp % array
- zx => grid % zx % array
- zz => grid % zz % array
- hx => grid % hx % array
- dss => grid % dss % array
-
- xCell => grid % xCell % array
- yCell => grid % yCell % array
-
- ppb => diag % pressure_base % array
- pb => diag % exner_base % array
- rb => diag % rho_base % array
- tb => diag % theta_base % array
- rtb => diag % rtheta_base % array
- p => diag % exner % array
- cqw => diag % cqw % array
-
- rho_zz => state % rho_zz % array
-
- pp => diag % pressure_p % array
- rr => diag % rho_p % array
- t => state % theta_m % array
- rt => diag % rtheta_p % array
- u => state % u % array
- ru => diag % ru % array
-
- scalars => state % scalars % array
-
- index_qv = state % index_qv
-
- scalars(:,:,:) = 0.
-
- call atm_initialize_advection_rk(grid)
- call atm_initialize_deformation_weights(grid)
-
- xnutr = 0.1
- zd = 10500.
-
- p0 = 1.e+05
- rcp = rgas/cp
- rcv = rgas/(cp-rgas)
-
- ! for hx computation
- xa = 5000. !SHP - should be changed based on grid distance
- xla = 4000.
- xc = maxval (grid % xCell % array(:))/2.
-
- ! metrics for hybrid coordinate and vertical stretching
- str = 1.0
- zt = 21000.
- dz = zt/float(nz1)
-! write(0,*) ' dz = ',dz
-
- do k=1,nz
-
-! sh(k) is the stretching specified for height surfaces
-
- zc(k) = zt*(real(k-1)*dz/zt)**str
-
-! to specify specific heights zc(k) for coordinate surfaces,
-! input zc(k)
-! zw(k) is the hieght of zeta surfaces
-! zw(k) = (k-1)*dz yields constant dzeta
-! and nonconstant dzeta/dz
-! zw(k) = sh(k)*zt yields nonconstant dzeta
-! and nearly constant dzeta/dz
-
-! zw(k) = float(k-1)*dz
- zw(k) = zc(k)
-!
-! ah(k) governs the transition between terrain-following
-! and pureheight coordinates
-! ah(k) = 0 is a terrain-following coordinate
-! ah(k) = 1 is a height coordinate
-
-! ah(k) = 1.-cos(.5*pii*(k-1)*dz/zt)**6
- ah(k) = 1.
-! write(0,*) ' k, zc, zw, ah ',k,zc(k),zw(k),ah(k)
- end do
- do k=1,nz1
- dzw (k) = zw(k+1)-zw(k)
- rdzw(k) = 1./dzw(k)
- zu(k ) = .5*(zw(k)+zw(k+1))
- end do
- do k=2,nz1
- dzu (k) = .5*(dzw(k)+dzw(k-1))
- rdzu(k) = 1./dzu(k)
- fzp (k) = .5* dzw(k )/dzu(k)
- fzm (k) = .5* dzw(k-1)/dzu(k)
- rdzwp(k) = dzw(k-1)/(dzw(k )*(dzw(k)+dzw(k-1)))
- rdzwm(k) = dzw(k )/(dzw(k-1)*(dzw(k)+dzw(k-1)))
- end do
-
-!********** how are we storing cf1, cf2 and cf3?
-
- d1 = .5*dzw(1)
- d2 = dzw(1)+.5*dzw(2)
- d3 = dzw(1)+dzw(2)+.5*dzw(3)
- !cf1 = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
- !cf2 = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
- !cf3 = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
-
- cof1 = (2.*dzu(2)+dzu(3))/(dzu(2)+dzu(3))*dzw(1)/dzu(2)
- cof2 = dzu(2) /(dzu(2)+dzu(3))*dzw(1)/dzu(3)
- cf1 = fzp(2) + cof1
- cf2 = fzm(2) - cof1 - cof2
- cf3 = cof2
-
- grid % cf1 % scalar = cf1
- grid % cf2 % scalar = cf2
- grid % cf3 % scalar = cf3
-
-! setting for terrain
- do iCell=1,grid % nCells
- xi = grid % xCell % array(iCell)
- !====1. for pure cosine mountain
- ! if(abs(xi-xc).ge.2.*xa) then
- ! hx(1,iCell) = 0.
- ! else
- ! hx(1,iCell) = hm*cos(.5*pii*(xi-xc)/(2.*xa))**2.
- ! end if
-
- !====2. for cosine mountain
- !if(abs(xi-xc).lt.xa) THEN
- ! hx(1,iCell) = hm*cos(pii*(xi-xc)/xla)**2. *cos(.5*pii*(xi-xc)/xa )**2.
- ! else
- ! hx(1,iCell) = 0.
- ! end if
-
- !====3. for shock mountain
- hx(1,iCell) = hm*exp(-((xi-xc)/xa)**2)*cos(pii*(xi-xc)/xla)**2.
-
- hx(nz,iCell) = zt
-
-!***** SHP -> get the temporary point information for the neighbor cell ->> should be changed!!!!!
- do i=1,grid % nCells
- !option 1
- !IF(yCell(i).eq.yCell(iCell).and.xCell(i).eq.xCell(iCell)-sqrt(3.)*grid % dcEdge % array(1)) next_cell(iCell,1) = i
- !IF(yCell(i).eq.yCell(iCell).and.xCell(i).eq.xCell(iCell)+sqrt(3.)*grid % dcEdge % array(1)) next_cell(iCell,2) = i
- !option 2
- next_cell(iCell,1) = iCell - 8 ! note ny=4
- next_cell(iCell,2) = iCell + 8 ! note ny=4
-
- if (xCell(iCell).le. 3.*grid % dcEdge % array(1)) then
- next_cell(iCell,1) = 1
- else if (xCell(iCell).ge. maxval(xCell(:))-3.*grid % dcEdge % array(1)) then
- next_cell(iCell,2) = 1
- end if
-
- end do
- enddo
-
- write(0,*) ' hx computation complete '
-
-
-! smoothing grid for the upper level >> but not propoer for parallel programing
- dzmin=.7
- do k=2,nz1
- sm = .25*min((zc(k)-zc(k-1))/dz,1.0_RKIND)
- do i=1,grid % nCells
- hx(k,i) = hx(k-1,i)
- end do
-
- do iter = 1,20 !iteration for smoothing
-
- do i=1,grid % nCells
- hxzt(i) = hx(k,i) + sm*(hx(k,next_cell(i,2))-2.*hx(k,i)+hx(k,next_cell(i,1)))
- end do
- dzh = zc(k) - zc(k-1)
- do i=1,grid % nCells
- dzht = zc(k)+hxzt(i) - zc(k-1)-hx(k-1,i)
- if(dzht.lt.dzh) dzh = dzht
- end do
-
- if(dzh.gt.dzmin*(zc(k)-zc(k-1))) then
- do i=1,grid % nCells
- hx(k,i) = hxzt(i)
- end do
- else
- goto 99 !SHP - this algorithm should be changed
- end if
-
- end do !end of iteration for smoothing
-99 write(0,*) "PASS-SHP"
- end do
-
- do iCell=1,grid % nCells
- do k=1,nz
- if (terrain_smooth) then
- zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(k,iCell)/zt)+hx(k,iCell)) &
- + (1.-ah(k)) * zc(k)
- else
- zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(1,iCell)/zt)+hx(1,iCell)) &
- + (1.-ah(k)) * zc(k)
- end if
- end do
- do k=1,nz1
- zz (k,iCell) = (zw(k+1)-zw(k))/(zgrid(k+1,iCell)-zgrid(k,iCell))
- end do
- end do
-
- do i=1, grid % nEdges
- iCell1 = grid % CellsOnEdge % array(1,i)
- iCell2 = grid % CellsOnEdge % array(2,i)
- do k=1,nz
- zx (k,i) = (zgrid(k,iCell2)-zgrid(k,iCell1)) / grid % dcEdge % array(i)
- end do
- end do
- do i=1, grid % nCells
- do k=1,nz1
- ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
- dss(k,i) = 0.
- ztemp = zgrid(k,i)
- if(ztemp.gt.zd+.1) then
- dss(k,i) = dss(k,i)+xnutr*sin(.5*pii*(ztemp-zd)/(zt-zd))**2
- end if
- end do
- enddo
-
- write(0,*) ' grid metrics setup complete '
-
-!
-! mountain wave initialization
-!
- !SHP-original
- !zinv = 1000.
- !SHP-schar case
- zinv = 3000.
-
- xn2 = 0.0001
- xn2m = 0.0000
- xn2l = 0.0001
-
- um = 10.
- us = 0.
-
- do i=1,grid % nCells
- do k=1,nz1
- ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
- tb(k,i) = t0*(1. + xn2m/gravity*ztemp)
- if(ztemp .le. zinv) then
- t (k,i) = t0*(1.+xn2l/gravity*ztemp)
- else
- t (k,i) = t0*(1.+xn2l/gravity*zinv+xn2/gravity*(ztemp-zinv))
- end if
- rh(k,i) = 0.
- end do
- end do
-
-! set the velocity field - we are on a plane here.
-
- do i=1, grid % nEdges
- cell1 = grid % CellsOnEdge % array(1,i)
- cell2 = grid % CellsOnEdge % array(2,i)
- if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,nz1
- ztemp = .25*( zgrid(k,cell1 )+zgrid(k+1,cell1 ) &
- +zgrid(k,cell2)+zgrid(k+1,cell2))
- u(k,i) = um
- if(i == 1 ) grid % u_init % array(k) = u(k,i) - us
-#ifdef ROTATED_GRID
- u(k,i) = sin(grid % angleEdge % array(i)) * (u(k,i) - us)
-#else
- u(k,i) = cos(grid % angleEdge % array(i)) * (u(k,i) - us)
-#endif
- end do
- end if
- end do
-
-!
-! reference sounding based on dry atmosphere
-!
- pitop = 1.-.5*dzw(1)*gravity/(cp*tb(1,1)*zz(1,1))
- do k=2,nz1
- pitop = pitop-dzu(k)*gravity/(cp*(fzm(k)*tb(k,1)+fzp(k)*tb(k-1,1)) &
- *(fzm(k)*zz(k,1)+fzp(k)*zz(k-1,1)))
- end do
- pitop = pitop-.5*dzw(nz1)*gravity/(cp*tb(nz1,1)*zz(nz1,1))
- ptopb = p0*pitop**(1./rcp)
-
- do i=1, grid % nCells
- pb(nz1,i) = pitop+.5*dzw(nz1)*gravity/(cp*tb(nz1,i)*zz(nz1,i))
- p (nz1,i) = pitop+.5*dzw(nz1)*gravity/(cp*t (nz1,i)*zz(nz1,i))
- do k=nz1-1,1,-1
- pb(k,i) = pb(k+1,i) + dzu(k+1)*gravity/(cp*.5*(tb(k,i)+tb(k+1,i)) &
- *.5*(zz(k,i)+zz(k+1,i)))
- p (k,i) = p (k+1,i) + dzu(k+1)*gravity/(cp*.5*(t (k,i)+t (k+1,i)) &
- *.5*(zz(k,i)+zz(k+1,i)))
- end do
- do k=1,nz1
- rb (k,i) = pb(k,i)**(1./rcv)/((rgas/p0)*tb(k,i)*zz(k,i))
- rtb(k,i) = rb(k,i)*tb(k,i)
- rr (k,i) = p (k,i)**(1./rcv)/((rgas/p0)*t (k,i)*zz(k,i))-rb(k,i)
- cqw(k,i) = 1.
- end do
- end do
-
- write(0,*) ' ***** base state sounding ***** '
- write(0,*) 'k pb p rb rtb rr tb t'
- do k=1,grid%nVertLevels
- write(0,'(i2,7(2x,f14.9))') k,pb(k,1),p(k,1),rb(k,1),rtb(k,1),rr(k,1),tb(k,1),t(k,1)
- end do
-
- scalars(index_qv,:,:) = 0.
-
-!-------------------------------------------------------------------
-! ITERATIONS TO CONVERGE MOIST SOUNDING
- do itr=1,30
- pitop = 1.-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1))
-
- do k=2,nz1
- pitop = pitop-dzu(k)*gravity/(cp*cqw(k,1)*(fzm(k)*t (k,1)+fzp(k)*t (k-1,1)) &
- *(fzm(k)*zz(k,1)+fzp(k)*zz(k-1,1)))
- end do
- pitop = pitop - .5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,1))/(cp*t(nz1,1)*zz(nz1,1))
- ptop = p0*pitop**(1./rcp)
-
- do i = 1, grid % nCells
-
- pp(nz1,i) = ptop-ptopb+.5*dzw(nz1)*gravity* &
- (rr(nz1,i)+(rr(nz1,i)+rb(nz1,i))*scalars(index_qv,nz1,i))
- do k=nz1-1,1,-1
- pp(k,i) = pp(k+1,i)+dzu(k+1)*gravity* &
- (fzm(k)*(rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i)) &
- +fzp(k)*(rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i)))
- end do
- do k=1,nz1
- rt(k,i) = (pp(k,i)/(rgas*zz(k,i)) &
- -rtb(k,i)*(p(k,i)-pb(k,i)))/p(k,i)
- p (k,i) = (zz(k,i)*(rgas/p0)*(rtb(k,i)+rt(k,i)))**rcv
- rr(k,i) = (rt(k,i)-rb(k,i)*(t(k,i)-tb(k,i)))/t(k,i)
- end do
-!
-! update water vapor mixing ratio from humitidty profile
-!
- do k=1,nz1
- temp = p(k,i)*t(k,i)
- pres = p0*p(k,i)**(1./rcp)
- qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
- scalars(index_qv,k,i) = min(0.014_RKIND,rh(k,i)*qvs)
- end do
-
- do k=1,nz1
- t (k,i) = t(k,i)*(1.+1.61*scalars(index_qv,k,i))
- end do
- do k=2,nz1
- cqw(k,i) = 1./(1.+.5*( scalars(index_qv,k-1,i) &
- +scalars(index_qv,k ,i)))
- end do
-
- end do ! loop over cells
-
- end do ! iteration loop
-!----------------------------------------------------------------------
-!
- write(0,*) ' *** sounding for the simulation ***'
- write(0,*) ' z theta pres qv rho_m u rr'
- do k=1,nz1
- write(6,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
- t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
- .01*p0*p(k,1)**(1./rcp), &
- 1000.*scalars(index_qv,k,1), &
- (rb(k,1)+rr(k,1))*(1.+scalars(index_qv,k,1)), &
- grid % u_init % array(k), rr(k,1)
- end do
-
- do i=1,grid % ncells
- do k=1,nz1
- rho_zz(k,i) = rb(k,i)+rr(k,i)
- end do
-
- do k=1,nz1
- grid % t_init % array(k,i) = t(k,i)
- end do
- end do
-
- do i=1,grid % nEdges
- cell1 = grid % CellsOnEdge % array(1,i)
- cell2 = grid % CellsOnEdge % array(2,i)
- if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,nz1
- ru (k,i) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*u(k,i)
- end do
- end if
- end do
-
-!
-! pre-calculation z-metric terms in omega eqn.
-!
- do iEdge = 1,grid % nEdges
- cell1 = CellsOnEdge(1,iEdge)
- cell2 = CellsOnEdge(2,iEdge)
- if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
-
- do k = 1, grid%nVertLevels
-
- if (config_theta_adv_order == 2) then
-
- z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
-
- else !theta_adv_order == 3 or 4
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
- do i=1, grid % nEdgesOnCell % array (cell1)
- if ( grid % CellsOnCell % array (i,cell1) > 0) &
- d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
- end do
- do i=1, grid % nEdgesOnCell % array (cell2)
- if ( grid % CellsOnCell % array (i,cell2) > 0) &
- d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
- - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- if (config_theta_adv_order == 3) then
- z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
- else
- z_edge3 = 0.
- end if
-
- end if
-
- zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
- zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
- zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
- zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
-
- end do
-
- end if
- end do
-
-! for including terrain
- state % w % array(:,:) = 0.0
- diag % rw % array(:,:) = 0.0
-
-!
-! calculation of omega, rw = zx * ru + zz * rw
-!
-
- do iEdge = 1,grid % nEdges
-
- cell1 = CellsOnEdge(1,iEdge)
- cell2 = CellsOnEdge(2,iEdge)
-
- if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
- do k = 2, grid%nVertLevels
- flux = (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))
- diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + (fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb(k,2,iEdge)*flux
- diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - (fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb(k,1,iEdge)*flux
-
- if (config_theta_adv_order ==3) then
- diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
- - sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order* &
- (fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb3(k,2,iEdge)*flux
- diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
- + sign(1.0_RKIND,ru(k,iEdge))*config_coef_3rd_order* &
- (fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb3(k,1,iEdge)*flux
- end if
-
- end do
- end if
-
- end do
-
- ! Compute w from rho_zz and rw
- do iCell=1,grid%nCells
- do k=2,grid%nVertLevels
- state % w % array(k,iCell) = diag % rw % array(k,iCell) &
- / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
- end do
- end do
-
-
- do iEdge=1,grid % nEdges
- grid % fEdge % array(iEdge) = 0.
- end do
-
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 0.
- end do
-
- !
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
- diag % v % array(:,:) = 0.0
- do iEdge = 1, grid%nEdges
- do i=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
- do k = 1, grid%nVertLevels
- diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
- end do
- end if
- end do
- end do
-
-! do k=1,grid%nVertLevels
-! write(0,*) ' k,u_init, t_init, qv_init ',k,grid % u_init % array(k),grid % t_init% array(k),grid % qv_init % array(k)
-! end do
-
- ! Compute rho and theta from rho_zz and theta_m
- do iCell=1,grid%nCells
- do k=1,grid%nVertLevels
- diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
- diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
- end do
- end do
-
- end subroutine atm_test_case_mtn_wave
-
-
-!----------------------------------------------------------------------------------------------------------
-
- real (kind=RKIND) function sphere_distance(lat1, lon1, lat2, lon2, radius)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
- ! sphere with given radius.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
- real (kind=RKIND) :: arg1
-
- arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
- cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
- sphere_distance = 2.*radius*asin(arg1)
-
- end function sphere_distance
-
-!--------------------------------------------------------------------
- real (kind=RKIND) function env_qv( z, temperature, pressure, rh_max )
-
- implicit none
- real (kind=RKIND) :: z, temperature, pressure, ztr, es, qvs, p0, rh_max
-
- p0 = 100000.
-
-! ztr = 5000.
-!
-! if(z .gt. ztr) then
-! env_qv = 0.
-! else
-! if(z.lt.2000.) then
-! env_qv = .5
-! else
-! env_qv = .5*(1.-(z-2000.)/(ztr-2000.))
-! end if
-! end if
-
- if (pressure .lt. 50000. ) then
- env_qv = 0.0
- else
- env_qv = (1.-((p0-pressure)/50000.)**1.25)
- end if
-
- env_qv = min(rh_max,env_qv)
-
-! env_qv is the relative humidity, turn it into mixing ratio
- if (temperature .gt. 273.15) then
- es = 1000.*0.6112*exp(17.67*(temperature-273.15)/(temperature-29.65))
- else
- es = 1000.*0.6112*exp(21.8745584*(temperature-273.16)/(temperature-7.66))
- end if
- qvs = (287.04/461.6)*es/(pressure-es)
-
- ! qvs = 380.*exp(17.27*(temperature-273.)/(temperature-36.))/pressure
-
- env_qv = env_qv*qvs
-
- end function env_qv
-
-end module atm_test_cases
Modified: branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_time_integration.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -1326,7 +1326,7 @@
wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
end do
end do
- do k=1,grid % nVertLevelsSolve
+ do k=1,grid % nVertLevels ! Could be nVertLevelsSolve?
do iScalar=1,s_old % num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
@@ -1355,7 +1355,7 @@
wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
enddo
- do k=1,grid % nVertLevelsSolve
+ do k=1,grid % nVertLevels ! Could be nVertLevelsSolve?
do iScalar=1,s_old % num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
@@ -1383,7 +1383,7 @@
wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
enddo
- do k=1,grid % nVertLevelsSolve
+ do k=1,grid % nVertLevels ! Could be nVertLevelsSolve?
do iScalar=1,s_old % num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
@@ -1865,12 +1865,11 @@
!SHP-curvature
logical, parameter :: curvature = .true.
- real (kind=RKIND), parameter :: omega_e = 7.29212e-05
real (kind=RKIND) :: r_earth
real (kind=RKIND), dimension(:,:), pointer :: ur_cell, vr_cell
real (kind=RKIND), parameter :: c_s = 0.125
-! real (kind=RKIND), parameter :: c_s = 0.25
+! real (kind=RKIND), parameter :: c_s = 0.25
real (kind=RKIND), dimension( grid % nVertLevels ) :: d_diag, d_off_diag, flux_arr
real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
logical :: delsq_horiz_mixing, newpx
@@ -1891,7 +1890,7 @@
!-----------
!SHP-curvature
- r_earth = a
+ r_earth = grid % sphere_radius
ur_cell => diag % uReconstructZonal % array
vr_cell => diag % uReconstructMeridional % array
@@ -2152,10 +2151,19 @@
tend_u(k,iEdge) = tend_u(k,iEdge) + rho_edge(k,iEdge)* (q(k) - (ke(k,cell2) - ke(k,cell1)) &
/ dcEdge(iEdge)) &
- u(k,iEdge)*0.5*(divergence_ru(k,cell1)+divergence_ru(k,cell2))
+ !SHP-curvature
+ if (curvature) then
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - 2.*omega_e*cos(grid % angleEdge % array(iEdge))*cos(grid % latEdge % array(iEdge)) &
- *.25*(rw(k,cell1)+rw(k+1,cell1)+rw(k,cell2)+rw(k+1,cell2)) &
- - u(k,iEdge)*.25*(rw(k+1,cell1)+rw(k,cell1)+rw(k,cell2)+rw(k+1,cell2))/r_earth
+ - 2.*omega*cos(grid % angleEdge % array(iEdge))*cos(grid % latEdge % array(iEdge)) &
+ *rho_edge(k,iEdge)*.25*(w(k,cell1)+w(k+1,cell1)+w(k,cell2)+w(k+1,cell2)) &
+ - u(k,iEdge)*.25*(w(k+1,cell1)+w(k,cell1)+w(k,cell2)+w(k+1,cell2)) &
+ *rho_edge(k,iEdge)/r_earth
+ !old-err.
+ !tend_u(k,iEdge) = tend_u(k,iEdge) &
+ ! - 2.*omega_e*cos(grid % angleEdge % array(iEdge))*cos(grid % latEdge % array(iEdge)) &
+ ! *.25*(rw(k,cell1)+rw(k+1,cell1)+rw(k,cell2)+rw(k+1,cell2)) &
+ ! - u(k,iEdge)*.25*(rw(k+1,cell1)+rw(k,cell1)+rw(k,cell2)+rw(k+1,cell2))/r_earth
+ end if
end do
end do
@@ -2506,12 +2514,19 @@
do iCell = 1, grid % nCellsSolve
do k=2,nVertLevels
- tend_w(k,iCell) = tend_w(k,iCell) &
- + rho_zz(k,iCell)*( (fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))**2. &
- +(fzm(k)*vr_cell(k,iCell)+fzp(k)*vr_cell(k-1,iCell))**2. )/r_earth &
- + 2.*omega_e*cos(grid % latCell % array(iCell))*rho_zz(k,iCell) &
- *(fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))
+ tend_w(k,iCell) = tend_w(k,iCell) + (rho_zz(k,iCell)*fzm(k)+rho_zz(k-1,iCell)*fzp(k))* &
+ ( (fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))**2. &
+ +(fzm(k)*vr_cell(k,iCell)+fzp(k)*vr_cell(k-1,iCell))**2. )/r_earth &
+ + 2.*omega*cos(grid % latCell % array(iCell)) &
+ *(fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell)) &
+ *(rho_zz(k,iCell)*fzm(k)+rho_zz(k-1,iCell)*fzp(k))
+ !old_err.
+ !tend_w(k,iCell) = tend_w(k,iCell) &
+ ! + rho_zz(k,iCell)*( (fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))**2. &
+ ! +(fzm(k)*vr_cell(k,iCell)+fzp(k)*vr_cell(k-1,iCell))**2. )/r_earth &
+ ! + 2.*omega_e*cos(grid % latCell % array(iCell))*rho_zz(k,iCell) &
+ ! *(fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))
end do
end do
Modified: branches/omp_blocks/openmp_test/src/framework/mpas_block_creator.F
===================================================================
--- branches/omp_blocks/openmp_test/src/framework/mpas_block_creator.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/framework/mpas_block_creator.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -1041,7 +1041,6 @@
block_ptr % mesh % nCellsArray(:) = nCellsCursor % array(:)
block_ptr % mesh % nEdgesArray(:) = nEdgesCursor % array(:)
block_ptr % mesh % nVerticesArray(:) = nVerticesCursor % array(:)
- block_ptr % mesh % nVertLevelsSolve = nVertLevels ! No vertical Decomposition yet...
! Set block's local id
block_ptr % localBlockID = localBlockID
Modified: branches/omp_blocks/openmp_test/src/framework/mpas_framework.F
===================================================================
--- branches/omp_blocks/openmp_test/src/framework/mpas_framework.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/framework/mpas_framework.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -13,18 +13,19 @@
contains
- subroutine mpas_framework_init(dminfo, domain)
+ subroutine mpas_framework_init(dminfo, domain, mpi_comm)
implicit none
type (dm_info), pointer :: dminfo
type (domain_type), pointer :: domain
+ integer, intent(in), optional :: mpi_comm
integer :: pio_num_iotasks
integer :: pio_stride
allocate(dminfo)
- call mpas_dmpar_init(dminfo)
+ call mpas_dmpar_init(dminfo, mpi_comm)
call mpas_read_namelist(dminfo)
Modified: branches/omp_blocks/openmp_test/src/framework/mpas_io_input.F
===================================================================
--- branches/omp_blocks/openmp_test/src/framework/mpas_io_input.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/framework/mpas_io_input.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -30,7 +30,6 @@
integer :: readCellStart, readCellEnd, nReadCells
integer :: readEdgeStart, readEdgeEnd, nReadEdges
integer :: readVertexStart, readVertexEnd, nReadVertices
- integer :: readVertLevelStart, readVertLevelEnd, nReadVertLevels
contains
@@ -147,10 +146,6 @@
call mpas_dmpar_get_index_range(domain % dminfo, 1, nVertices, readVertexStart, readVertexEnd)
nReadVertices = readVertexEnd - readVertexStart + 1
- readVertLevelStart = 1
- readVertLevelEnd = nVertLevels
- nReadVertLevels = nVertLevels
-
allocate(readingBlock)
readingBlock % domain => domain
readingBlock % blockID = domain % dminfo % my_proc_id
@@ -245,7 +240,6 @@
do while (associated(block_ptr))
block_ptr % mesh % sphere_radius = domain % blocklist % mesh % sphere_radius
block_ptr % mesh % on_a_sphere = domain % blocklist % mesh % on_a_sphere
- block_ptr % mesh % nVertLevelsSolve = domain % blocklist % mesh % nVertLevelsSolve ! No vertical decomp yet...
! Link the sendList and recvList pointers in each field type to the appropriate lists
! in parinfo, e.g., cellsToSend and cellsToRecv; in future, it can also be extended to
@@ -279,7 +273,7 @@
! process
! 2) All processes then send the global indices that were read to the
! processes that own those indices based on
- ! {send,recv}{Cell,Edge,Vertex,VertLevel}List
+ ! {send,recv}{Cell,Edge,Vertex}List
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call mpas_read_and_distribute_fields(input_obj)
Modified: branches/omp_blocks/openmp_test/src/framework/mpas_timer.F
===================================================================
--- branches/omp_blocks/openmp_test/src/framework/mpas_timer.F 2012-09-06 21:16:11 UTC (rev 2144)
+++ branches/omp_blocks/openmp_test/src/framework/mpas_timer.F 2012-09-07 15:07:43 UTC (rev 2145)
@@ -46,6 +46,10 @@
integer :: clock, hz, usecs
+#ifdef MPAS_TAU
+ call tau_start(timer_name)
+#endif
+
timer_added = .false.
timer_found = .false.
@@ -170,6 +174,10 @@
real (kind=RKIND) :: time_temp
logical :: timer_found, string_equal, check_flag
integer :: clock, hz, usecs
+
+#ifdef MPAS_TAU
+ call tau_stop(timer_name)
+#endif
timer_found = .false.
</font>
</pre>