<p><b>duda</b> 2010-11-03 12:27:30 -0600 (Wed, 03 Nov 2010)</p><p>BRANCH COMMIT<br>
<br>
Add new initialization core, init_nhyd_atmos, for non-hydrostatic model.<br>
This core only generates initial conditions, which can be used by the<br>
non-hydrostatic model by specifying test case 0 (use ICs from input file).<br>
<br>
The initialization core supports all of the test cases available in<br>
the non-hydrostatic model.<br>
<br>
<br>
A src/core_init_nhyd_atmos<br>
A src/core_init_nhyd_atmos/module_core.F<br>
A src/core_init_nhyd_atmos/module_advection.F<br>
A src/core_init_nhyd_atmos/module_test_cases.F<br>
A src/core_init_nhyd_atmos/Registry<br>
A src/core_init_nhyd_atmos/Makefile<br>
</p><hr noshade><pre><font color="gray">Added: branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/Makefile
===================================================================
--- branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/Makefile (rev 0)
+++ branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/Makefile 2010-11-03 18:27:30 UTC (rev 592)
@@ -0,0 +1,27 @@
+.SUFFIXES: .F .o
+
+OBJS = module_core.o \
+ module_test_cases.o \
+ module_advection.o
+
+all: core_hyd
+
+core_hyd: $(OBJS)
+ ar -ru libdycore.a $(OBJS)
+
+module_test_cases.o: module_advection.o
+
+module_advection.o:
+
+module_core.o: module_advection.o module_test_cases.o
+
+clean:
+ $(RM) *.o *.mod *.f90 libdycore.a
+
+.F.o:
+ $(RM) $@ $*.mod
+ $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
+ $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
+
+.c.o:
+ $(CC) $(CFLAGS) $(CPPFLAGS) $(CPPINCLUDES) -c $<
Added: branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/Registry (rev 0)
+++ branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/Registry 2010-11-03 18:27:30 UTC (rev 592)
@@ -0,0 +1,188 @@
+#
+# namelist type namelist_record name default_value
+#
+namelist integer nhyd_model config_test_case 5
+namelist character nhyd_model config_time_integration SRK3
+namelist real nhyd_model config_dt 172.8
+namelist integer nhyd_model config_ntimesteps 7500
+namelist integer nhyd_model config_output_interval 500
+namelist character nhyd_model config_horiz_mixing 2d_smagorinsky
+namelist real nhyd_model config_h_mom_eddy_visc2 0.0
+namelist real nhyd_model config_h_mom_eddy_visc4 0.0
+namelist real nhyd_model config_v_mom_eddy_visc2 0.0
+namelist real nhyd_model config_h_theta_eddy_visc2 0.0
+namelist real nhyd_model config_h_theta_eddy_visc4 0.0
+namelist real nhyd_model config_v_theta_eddy_visc2 0.0
+namelist integer nhyd_model config_number_of_sub_steps 4
+namelist integer nhyd_model config_w_adv_order 2
+namelist integer nhyd_model config_theta_adv_order 2
+namelist integer nhyd_model config_scalar_adv_order 2
+namelist integer nhyd_model config_u_vadv_order 2
+namelist integer nhyd_model config_w_vadv_order 2
+namelist integer nhyd_model config_theta_vadv_order 2
+namelist integer nhyd_model config_scalar_vadv_order 2
+namelist real nhyd_model config_coef_3rd_order 1.0
+namelist logical nhyd_model config_scalar_advection true
+namelist logical nhyd_model config_positive_definite false
+namelist logical nhyd_model config_monotonic true
+namelist logical nhyd_model config_mix_full true
+namelist real nhyd_model config_len_disp 0.
+namelist integer nhyd_model config_mp_physics 0.
+namelist real nhyd_model config_epssm 0.1
+namelist real nhyd_model config_smdiv 0.1
+namelist integer dimensions config_nvertlevels 26
+namelist character io config_input_name grid.nc
+namelist character io config_output_name output.nc
+namelist character io config_restart_name restart.nc
+namelist integer restart config_restart_interval 0
+namelist logical restart config_do_restart false
+namelist real restart config_restart_time 172800.0
+
+#
+# dim type name_in_file name_in_code
+#
+dim nCells nCells
+dim nEdges nEdges
+dim maxEdges maxEdges
+dim maxEdges2 maxEdges2
+dim nVertices nVertices
+dim TWO 2
+dim THREE 3
+dim vertexDegree vertexDegree
+dim FIFTEEN 15
+dim TWENTYONE 21
+dim R3 3
+dim nVertLevels namelist:config_nvertlevels
+dim nVertLevelsP1 nVertLevels+1
+
+#
+# var type name_in_file ( dims ) iro- name_in_code super-array array_class
+#
+var persistent real xtime ( Time ) 2 o xtime state - -
+
+# horizontal grid structure
+
+var persistent real latCell ( nCells ) 0 io latCell mesh - -
+var persistent real lonCell ( nCells ) 0 io lonCell mesh - -
+var persistent real xCell ( nCells ) 0 io xCell mesh - -
+var persistent real yCell ( nCells ) 0 io yCell mesh - -
+var persistent real zCell ( nCells ) 0 io zCell mesh - -
+var persistent integer indexToCellID ( nCells ) 0 io indexToCellID mesh - -
+
+var persistent real latEdge ( nEdges ) 0 io latEdge mesh - -
+var persistent real lonEdge ( nEdges ) 0 io lonEdge mesh - -
+var persistent real xEdge ( nEdges ) 0 io xEdge mesh - -
+var persistent real yEdge ( nEdges ) 0 io yEdge mesh - -
+var persistent real zEdge ( nEdges ) 0 io zEdge mesh - -
+var persistent integer indexToEdgeID ( nEdges ) 0 io indexToEdgeID mesh - -
+
+var persistent real latVertex ( nVertices ) 0 io latVertex mesh - -
+var persistent real lonVertex ( nVertices ) 0 io lonVertex mesh - -
+var persistent real xVertex ( nVertices ) 0 io xVertex mesh - -
+var persistent real yVertex ( nVertices ) 0 io yVertex mesh - -
+var persistent real zVertex ( nVertices ) 0 io zVertex mesh - -
+var persistent integer indexToVertexID ( nVertices ) 0 io indexToVertexID mesh - -
+
+var persistent integer cellsOnEdge ( TWO nEdges ) 0 io cellsOnEdge mesh - -
+var persistent integer nEdgesOnCell ( nCells ) 0 io nEdgesOnCell mesh - -
+var persistent integer nEdgesOnEdge ( nEdges ) 0 io nEdgesOnEdge mesh - -
+var persistent integer edgesOnCell ( maxEdges nCells ) 0 io edgesOnCell mesh - -
+var persistent integer edgesOnEdge ( maxEdges2 nEdges ) 0 io edgesOnEdge mesh - -
+
+var persistent real weightsOnEdge ( maxEdges2 nEdges ) 0 io weightsOnEdge mesh - -
+var persistent real dvEdge ( nEdges ) 0 io dvEdge mesh - -
+var persistent real dcEdge ( nEdges ) 0 io dcEdge mesh - -
+var persistent real angleEdge ( nEdges ) 0 io angleEdge mesh - -
+var persistent real areaCell ( nCells ) 0 io areaCell mesh - -
+var persistent real areaTriangle ( nVertices ) 0 io areaTriangle mesh - -
+
+var persistent real edgeNormalVectors ( R3 nEdges ) 0 io edgeNormalVectors mesh - -
+var persistent real localVerticalUnitVectors ( R3 nCells ) 0 io localVerticalUnitVectors mesh - -
+var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 io cellTangentPlane mesh - -
+
+var persistent integer cellsOnCell ( maxEdges nCells ) 0 io cellsOnCell mesh - -
+var persistent integer verticesOnCell ( maxEdges nCells ) 0 io verticesOnCell mesh - -
+var persistent integer verticesOnEdge ( TWO nEdges ) 0 io verticesOnEdge mesh - -
+var persistent integer edgesOnVertex ( vertexDegree nVertices ) 0 io edgesOnVertex mesh - -
+var persistent integer cellsOnVertex ( vertexDegree nVertices ) 0 io cellsOnVertex mesh - -
+var persistent real kiteAreasOnVertex ( vertexDegree nVertices ) 0 io kiteAreasOnVertex mesh - -
+var persistent real fEdge ( nEdges ) 0 io fEdge mesh - -
+var persistent real fVertex ( nVertices ) 0 io fVertex mesh - -
+
+# some solver scalar coefficients
+
+# coefficients for vertical extrapolation to the surface
+var persistent real cf1 ( ) 0 io cf1 mesh - -
+var persistent real cf2 ( ) 0 io cf2 mesh - -
+var persistent real cf3 ( ) 0 io cf3 mesh - -
+
+# description of the vertical grid structure
+
+var persistent real hx ( nVertLevelsP1 nCells ) 0 io hx mesh - -
+var persistent real zgrid ( nVertLevelsP1 nCells ) 0 io zgrid mesh - -
+var persistent real rdzw ( nVertLevels ) 0 io rdzw mesh - -
+var persistent real dzu ( nVertLevels ) 0 io dzu mesh - -
+var persistent real rdzu ( nVertLevels ) 0 io rdzu mesh - -
+var persistent real fzm ( nVertLevels ) 0 io fzm mesh - -
+var persistent real fzp ( nVertLevels ) 0 io fzp mesh - -
+var persistent real zx ( nVertLevelsP1 nEdges ) 0 io zx mesh - -
+var persistent real zz ( nVertLevelsP1 nCells ) 0 io zz mesh - -
+var persistent real zf ( nVertLevelsP1 TWO nEdges ) 0 io zf mesh - -
+var persistent real zf3 ( nVertLevelsP1 TWO nEdges ) 0 io zf3 mesh - -
+var persistent real zb ( nVertLevelsP1 TWO nEdges ) 0 io zb mesh - -
+var persistent real zb3 ( nVertLevelsP1 TWO nEdges ) 0 io zb3 mesh - -
+
+# W-Rayleigh-damping coefficient
+
+var persistent real dss ( nVertLevels nCells ) 0 io dss mesh - -
+
+# Prognostic variables: read from input, saved in restart, and written to output
+var persistent real u ( nVertLevels nEdges Time ) 2 io u state - -
+var persistent real w ( nVertLevelsP1 nCells Time ) 2 io w state - -
+var persistent real rho ( nVertLevels nCells Time ) 2 io rho state - -
+var persistent real theta ( nVertLevels nCells Time ) 2 io theta state - -
+var persistent real qv ( nVertLevels nCells Time ) 2 io qv state scalars moist
+var persistent real qc ( nVertLevels nCells Time ) 2 io qc state scalars moist
+var persistent real qr ( nVertLevels nCells Time ) 2 io qr state scalars moist
+
+# state variables diagnosed from prognostic state
+var persistent real pressure_p ( nVertLevels nCells Time ) 1 - pressure_p diag - -
+
+var persistent real u_init ( nVertLevels ) 0 io u_init mesh - -
+var persistent real t_init ( nVertLevels nCells ) 0 io t_init mesh - -
+var persistent real qv_init ( nVertLevels ) 0 io qv_init mesh - -
+
+# Diagnostic fields: only written to output
+var persistent real v ( nVertLevels nEdges Time ) 1 o v diag - -
+var persistent real uReconstructX ( nVertLevels nCells Time ) 1 o uReconstructX diag - -
+var persistent real uReconstructY ( nVertLevels nCells Time ) 1 o uReconstructY diag - -
+var persistent real uReconstructZ ( nVertLevels nCells Time ) 1 o uReconstructZ diag - -
+var persistent real uReconstructZonal ( nVertLevels nCells Time ) 1 o uReconstructZonal diag - -
+var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 1 o uReconstructMeridional diag - -
+
+var persistent real exner ( nVertLevels nCells Time ) 1 - exner diag - -
+var persistent real exner_base ( nVertLevels nCells Time ) 1 io exner_base diag - -
+var persistent real rtheta_base ( nVertLevels nCells Time ) 1 - rtheta_base diag - -
+var persistent real pressure_base ( nVertLevels nCells Time ) 1 io pressure_base diag - -
+var persistent real rho_base ( nVertLevels nCells Time ) 1 io rho_base diag - -
+var persistent real theta_base ( nVertLevels nCells Time ) 1 io theta_base diag - -
+
+var persistent real cqw ( nVertLevels nCells Time ) 1 - cqw diag - -
+
+# coupled variables needed by the solver, but not output...
+var persistent real ru ( nVertLevels nEdges Time ) 1 - ru diag - -
+var persistent real rw ( nVertLevelsP1 nCells Time ) 1 - rw diag - -
+var persistent real rtheta_p ( nVertLevels nCells Time ) 1 - rtheta_p diag - -
+var persistent real rho_p ( nVertLevels nCells Time ) 1 - rho_p diag - -
+
+# Space needed for advection
+var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 io deriv_two mesh - -
+var persistent integer advCells ( TWENTYONE nCells ) 0 io advCells mesh - -
+
+# Space needed for deformation calculation weights
+var persistent real defc_a ( maxEdges nCells ) 0 io defc_a mesh - -
+var persistent real defc_b ( maxEdges nCells ) 0 io defc_b mesh - -
+
+# Arrays required for reconstruction of velocity field
+var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 io coeffs_reconstruct mesh - -
+
Added: branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_advection.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_advection.F (rev 0)
+++ branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_advection.F 2010-11-03 18:27:30 UTC (rev 592)
@@ -0,0 +1,933 @@
+module advection
+
+ use grid_types
+ use configure
+ use constants
+
+
+ contains
+
+
+ subroutine initialize_advection_rk( grid )
+
+!
+! compute the cell coefficients for the polynomial fit.
+! this is performed during setup for model integration.
+! WCS, 31 August 2009
+!
+ implicit none
+
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ integer, dimension(:,:), pointer :: advCells
+
+! local variables
+
+ real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
+ real (kind=RKIND), dimension(grid % nEdges) :: xe, ye
+ real (kind=RKIND), dimension(grid % nCells) :: theta_abs
+
+ real (kind=RKIND), dimension(25) :: xc, yc, zc ! cell center coordinates
+ real (kind=RKIND), dimension(25) :: thetav, thetat, dl_sphere
+ real (kind=RKIND) :: xm, ym, zm, dl, xec, yec, zec
+ real (kind=RKIND) :: thetae_tmp, xe_tmp, ye_tmp
+ real (kind=RKIND) :: xv1, xv2, yv1, yv2, zv1, zv2
+ integer :: i, j, k, ip1, ip2, m, n, ip1a, ii
+ integer :: iCell, iEdge
+ real (kind=RKIND) :: pii
+ real (kind=RKIND) :: x0, y0, x1, y1, x2, y2, x3, y3, x4, y4, x5, y5
+ real (kind=RKIND) :: pdx1, pdx2, pdx3, pdy1, pdy2, pdy3, dx1, dx2, dy1, dy2
+ real (kind=RKIND) :: angv1, angv2, dl1, dl2
+ real (kind=RKIND), dimension(25) :: dxe, dye, x2v, y2v, xp, yp
+
+ real (kind=RKIND) :: amatrix(25,25), bmatrix(25,25), wmatrix(25,25)
+ real (kind=RKIND) :: length_scale
+ integer :: ma,na, cell_add, mw, nn
+ integer, dimension(25) :: cell_list
+
+
+ integer :: cell1, cell2
+ integer, parameter :: polynomial_order = 2
+! logical, parameter :: debug = .true.
+ logical, parameter :: debug = .false.
+! logical, parameter :: least_squares = .false.
+ logical, parameter :: least_squares = .true.
+ logical :: add_the_cell, do_the_cell
+
+ logical, parameter :: reset_poly = .true.
+
+ real (kind=RKIND) :: rcell, cos2t, costsint, sin2t
+ real (kind=RKIND), dimension(grid%maxEdges) :: angle_2d
+
+!---
+
+ pii = 2.*asin(1.0)
+
+ advCells => grid % advCells % array
+ deriv_two => grid % deriv_two % array
+ deriv_two(:,:,:) = 0.
+
+ do iCell = 1, grid % nCells ! is this correct? - we need first halo cell also...
+
+ cell_list(1) = iCell
+ do i=2, grid % nEdgesOnCell % array(iCell)+1
+ cell_list(i) = grid % CellsOnCell % array(i-1,iCell)
+ end do
+ n = grid % nEdgesOnCell % array(iCell) + 1
+
+ if ( polynomial_order > 2 ) then
+ do i=2,grid % nEdgesOnCell % array(iCell) + 1
+ do j=1,grid % nEdgesOnCell % array ( cell_list(i) )
+ cell_add = grid % CellsOnCell % array (j,cell_list(i))
+ add_the_cell = .true.
+ do k=1,n
+ if ( cell_add == cell_list(k) ) add_the_cell = .false.
+ end do
+ if (add_the_cell) then
+ n = n+1
+ cell_list(n) = cell_add
+ end if
+ end do
+ end do
+ end if
+
+ advCells(1,iCell) = n
+
+! check to see if we are reaching outside the halo
+
+ do_the_cell = .true.
+ do i=1,n
+ if (cell_list(i) > grid % nCells) do_the_cell = .false.
+ end do
+
+
+ if ( .not. do_the_cell ) cycle
+
+
+! compute poynomial fit for this cell if all needed neighbors exist
+ if ( grid % on_a_sphere ) then
+
+ 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
+ end do
+
+ theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
+ xc(2), yc(2), zc(2), &
+ 0., 0., 1. )
+
+! angles from cell center to neighbor centers (thetav)
+
+ do i=1,n-1
+
+ ip2 = i+2
+ if (ip2 > n) ip2 = 2
+
+ thetav(i) = sphere_angle( xc(1), yc(1), zc(1), &
+ 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) )
+ end do
+
+ length_scale = 1.
+ do i=1,n-1
+ dl_sphere(i) = dl_sphere(i)/length_scale
+ end do
+
+! thetat(1) = 0. ! this defines the x direction, cell center 1 ->
+ thetat(1) = theta_abs(iCell) ! this defines the x direction, longitude line
+ do i=2,n-1
+ thetat(i) = thetat(i-1) + thetav(i-1)
+ end do
+
+ do i=1,n-1
+ xp(i) = cos(thetat(i)) * dl_sphere(i)
+ yp(i) = sin(thetat(i)) * dl_sphere(i)
+ end do
+
+ else ! On an x-y plane
+
+ do i=1,n-1
+
+ angle_2d(i) = grid%angleEdge%array(grid % EdgesOnCell % array(i,iCell))
+ iEdge = grid % EdgesOnCell % array(i,iCell)
+ if ( iCell .ne. grid % CellsOnEdge % array(1,iEdge)) &
+ angle_2d(i) = angle_2d(i) - pii
+
+! xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
+! yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+
+ xp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * cos(angle_2d(i))
+ yp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * sin(angle_2d(i))
+
+ end do
+
+ end if
+
+
+ ma = n-1
+ mw = grid % nEdgesOnCell % array (iCell)
+
+ bmatrix = 0.
+ amatrix = 0.
+ wmatrix = 0.
+
+ if (polynomial_order == 2) then
+ na = 6
+ ma = ma+1
+
+ amatrix(1,1) = 1.
+ wmatrix(1,1) = 1.
+ do i=2,ma
+ amatrix(i,1) = 1.
+ amatrix(i,2) = xp(i-1)
+ amatrix(i,3) = yp(i-1)
+ amatrix(i,4) = xp(i-1)**2
+ amatrix(i,5) = xp(i-1) * yp(i-1)
+ amatrix(i,6) = yp(i-1)**2
+
+ wmatrix(i,i) = 1.
+ end do
+
+ else if (polynomial_order == 3) then
+ na = 10
+ ma = ma+1
+
+ amatrix(1,1) = 1.
+ wmatrix(1,1) = 1.
+ do i=2,ma
+ amatrix(i,1) = 1.
+ amatrix(i,2) = xp(i-1)
+ amatrix(i,3) = yp(i-1)
+
+ amatrix(i,4) = xp(i-1)**2
+ amatrix(i,5) = xp(i-1) * yp(i-1)
+ amatrix(i,6) = yp(i-1)**2
+
+ amatrix(i,7) = xp(i-1)**3
+ amatrix(i,8) = yp(i-1) * (xp(i-1)**2)
+ amatrix(i,9) = xp(i-1) * (yp(i-1)**2)
+ amatrix(i,10) = yp(i-1)**3
+
+ wmatrix(i,i) = 1.
+
+ end do
+
+ else
+ na = 15
+ ma = ma+1
+
+ amatrix(1,1) = 1.
+ wmatrix(1,1) = 1.
+ do i=2,ma
+ amatrix(i,1) = 1.
+ amatrix(i,2) = xp(i-1)
+ amatrix(i,3) = yp(i-1)
+
+ amatrix(i,4) = xp(i-1)**2
+ amatrix(i,5) = xp(i-1) * yp(i-1)
+ amatrix(i,6) = yp(i-1)**2
+
+ amatrix(i,7) = xp(i-1)**3
+ amatrix(i,8) = yp(i-1) * (xp(i-1)**2)
+ amatrix(i,9) = xp(i-1) * (yp(i-1)**2)
+ amatrix(i,10) = yp(i-1)**3
+
+ amatrix(i,11) = xp(i-1)**4
+ amatrix(i,12) = yp(i-1) * (xp(i-1)**3)
+ amatrix(i,13) = (xp(i-1)**2)*(yp(i-1)**2)
+ amatrix(i,14) = xp(i-1) * (yp(i-1)**3)
+ amatrix(i,15) = yp(i-1)**4
+
+ wmatrix(i,i) = 1.
+
+ end do
+
+ do i=1,mw
+ wmatrix(i,i) = 1.
+ end do
+
+ end if
+
+ call poly_fit_2( amatrix, bmatrix, wmatrix, ma, na, 25 )
+
+ do i=1,grid % nEdgesOnCell % array (iCell)
+ ip1 = i+1
+ 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
+
+ if ( grid % on_a_sphere ) then
+ call arc_bisect( xv1, yv1, zv1, &
+ xv2, yv2, zv2, &
+ xec, yec, zec )
+
+ thetae_tmp = sphere_angle( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1), &
+ xec, yec, zec )
+ thetae_tmp = thetae_tmp + thetat(i)
+ if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
+ thetae(1,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ else
+ thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ end if
+! else
+!
+! xe(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (xv1 + xv2)
+! ye(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (yv1 + yv2)
+
+ end if
+
+ end do
+
+! fill second derivative stencil for rk advection
+
+ do i=1, grid % nEdgesOnCell % array (iCell)
+ iEdge = grid % EdgesOnCell % array (i,iCell)
+
+
+ if ( grid % on_a_sphere ) then
+ if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
+
+ cos2t = cos(thetae(1,grid % EdgesOnCell % array (i,iCell)))
+ sin2t = sin(thetae(1,grid % EdgesOnCell % array (i,iCell)))
+ costsint = cos2t*sin2t
+ cos2t = cos2t**2
+ sin2t = sin2t**2
+
+ do j=1,n
+ deriv_two(j,1,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ + 2.*costsint*bmatrix(5,j) &
+ + 2.*sin2t*bmatrix(6,j)
+ end do
+ else
+
+ cos2t = cos(thetae(2,grid % EdgesOnCell % array (i,iCell)))
+ sin2t = sin(thetae(2,grid % EdgesOnCell % array (i,iCell)))
+ costsint = cos2t*sin2t
+ cos2t = cos2t**2
+ sin2t = sin2t**2
+
+ do j=1,n
+ deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ + 2.*costsint*bmatrix(5,j) &
+ + 2.*sin2t*bmatrix(6,j)
+ end do
+ end if
+
+ else
+
+ cos2t = cos(angle_2d(i))
+ sin2t = sin(angle_2d(i))
+ costsint = cos2t*sin2t
+ cos2t = cos2t**2
+ sin2t = sin2t**2
+
+! do j=1,n
+!
+! deriv_two(j,1,iEdge) = 2.*xe(iEdge)*xe(iEdge)*bmatrix(4,j) &
+! + 2.*xe(iEdge)*ye(iEdge)*bmatrix(5,j) &
+! + 2.*ye(iEdge)*ye(iEdge)*bmatrix(6,j)
+! end do
+
+ if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
+ do j=1,n
+ deriv_two(j,1,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ + 2.*costsint*bmatrix(5,j) &
+ + 2.*sin2t*bmatrix(6,j)
+ end do
+ else
+ do j=1,n
+ deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ + 2.*costsint*bmatrix(5,j) &
+ + 2.*sin2t*bmatrix(6,j)
+ end do
+ end if
+
+ end if
+ end do
+
+ end do ! end of loop over cells
+
+ if (debug) stop
+
+
+! write(0,*) ' check for deriv2 coefficients, iEdge 4 '
+!
+! iEdge = 4
+! j = 1
+! iCell = grid % cellsOnEdge % array(1,iEdge)
+! write(0,*) ' j, icell, coef ',j,iCell,deriv_two(j,1,iEdge)
+! do j=2,7
+! write(0,*) ' j, icell, coef ',j,grid % CellsOnCell % array(j-1,iCell),deriv_two(j,1,iEdge)
+! end do
+!
+! j = 1
+! iCell = grid % cellsOnEdge % array(2,iEdge)
+! write(0,*) ' j, icell, coef ',j,iCell,deriv_two(j,2,iEdge)
+! do j=2,7
+! write(0,*) ' j, icell, coef ',j,grid % CellsOnCell % array(j-1,iCell),deriv_two(j,2,iEdge)
+! end do
+! stop
+
+ end subroutine initialize_advection_rk
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! FUNCTION SPHERE_ANGLE
+ !
+ ! Computes the angle between arcs AB and AC, given points A, B, and C
+ ! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz
+
+ real (kind=RKIND) :: a, b, c ! Side lengths of spherical triangle ABC
+
+ real (kind=RKIND) :: ABx, ABy, ABz ! The components of the vector AB
+ real (kind=RKIND) :: mAB ! The magnitude of AB
+ real (kind=RKIND) :: ACx, ACy, ACz ! The components of the vector AC
+ real (kind=RKIND) :: mAC ! The magnitude of AC
+
+ real (kind=RKIND) :: Dx ! The i-components of the cross product AB x AC
+ real (kind=RKIND) :: Dy ! The j-components of the cross product AB x AC
+ real (kind=RKIND) :: Dz ! The k-components of the cross product AB x AC
+
+ real (kind=RKIND) :: s ! Semiperimeter of the triangle
+ real (kind=RKIND) :: sin_angle
+
+ a = acos(max(min(bx*cx + by*cy + bz*cz,1.0),-1.0)) ! Eqn. (3)
+ b = acos(max(min(ax*cx + ay*cy + az*cz,1.0),-1.0)) ! Eqn. (2)
+ c = acos(max(min(ax*bx + ay*by + az*bz,1.0),-1.0)) ! Eqn. (1)
+
+ ABx = bx - ax
+ ABy = by - ay
+ ABz = bz - az
+
+ ACx = cx - ax
+ ACy = cy - ay
+ ACz = cz - az
+
+ Dx = (ABy * ACz) - (ABz * ACy)
+ Dy = -((ABx * ACz) - (ABz * ACx))
+ Dz = (ABx * ACy) - (ABy * ACx)
+
+ s = 0.5*(a + b + c)
+! sin_angle = sqrt((sin(s-b)*sin(s-c))/(sin(b)*sin(c))) ! Eqn. (28)
+ sin_angle = sqrt(min(1.,max(0.,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
+
+ if ((Dx*ax + Dy*ay + Dz*az) >= 0.0) then
+ sphere_angle = 2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ else
+ sphere_angle = -2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ end if
+
+ end function sphere_angle
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! FUNCTION PLANE_ANGLE
+ !
+ ! Computes the angle between vectors AB and AC, given points A, B, and C, and
+ ! a vector (u,v,w) normal to the plane.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w
+
+ real (kind=RKIND) :: ABx, ABy, ABz ! The components of the vector AB
+ real (kind=RKIND) :: mAB ! The magnitude of AB
+ real (kind=RKIND) :: ACx, ACy, ACz ! The components of the vector AC
+ real (kind=RKIND) :: mAC ! The magnitude of AC
+
+ real (kind=RKIND) :: Dx ! The i-components of the cross product AB x AC
+ real (kind=RKIND) :: Dy ! The j-components of the cross product AB x AC
+ real (kind=RKIND) :: Dz ! The k-components of the cross product AB x AC
+
+ real (kind=RKIND) :: cos_angle
+
+ ABx = bx - ax
+ ABy = by - ay
+ ABz = bz - az
+ mAB = sqrt(ABx**2.0 + ABy**2.0 + ABz**2.0)
+
+ ACx = cx - ax
+ ACy = cy - ay
+ ACz = cz - az
+ mAC = sqrt(ACx**2.0 + ACy**2.0 + ACz**2.0)
+
+
+ Dx = (ABy * ACz) - (ABz * ACy)
+ Dy = -((ABx * ACz) - (ABz * ACx))
+ Dz = (ABx * ACy) - (ABy * ACx)
+
+ cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
+
+ if ((Dx*u + Dy*v + Dz*w) >= 0.0) then
+ plane_angle = acos(max(min(cos_angle,1.0),-1.0))
+ else
+ plane_angle = -acos(max(min(cos_angle,1.0),-1.0))
+ end if
+
+ end function plane_angle
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! FUNCTION ARC_LENGTH
+ !
+ ! Returns the length of the great circle arc from A=(ax, ay, az) to
+ ! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
+ ! same sphere centered at the origin.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real function arc_length(ax, ay, az, bx, by, bz)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz
+
+ real (kind=RKIND) :: r, c
+ real (kind=RKIND) :: cx, cy, cz
+
+ cx = bx - ax
+ cy = by - ay
+ cz = bz - az
+
+! r = ax*ax + ay*ay + az*az
+! c = cx*cx + cy*cy + cz*cz
+!
+! arc_length = sqrt(r) * acos(1.0 - c/(2.0*r))
+
+ r = sqrt(ax*ax + ay*ay + az*az)
+ c = sqrt(cx*cx + cy*cy + cz*cz)
+! arc_length = sqrt(r) * 2.0 * asin(c/(2.0*r))
+ arc_length = r * 2.0 * asin(c/(2.0*r))
+
+ end function arc_length
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE ARC_BISECT
+ !
+ ! Returns the point C=(cx, cy, cz) that bisects the great circle arc from
+ ! A=(ax, ay, az) to B=(bx, by, bz). It is assumed that A and B lie on the
+ ! surface of a sphere centered at the origin.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine arc_bisect(ax, ay, az, bx, by, bz, cx, cy, cz)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz
+ real (kind=RKIND), intent(out) :: cx, cy, cz
+
+ real (kind=RKIND) :: r ! Radius of the sphere
+ real (kind=RKIND) :: d
+
+ r = sqrt(ax*ax + ay*ay + az*az)
+
+ cx = 0.5*(ax + bx)
+ cy = 0.5*(ay + by)
+ cz = 0.5*(az + bz)
+
+ if (cx == 0. .and. cy == 0. .and. cz == 0.) then
+ write(0,*) 'Error: arc_bisect: A and B are diametrically opposite'
+ else
+ d = sqrt(cx*cx + cy*cy + cz*cz)
+ cx = r * cx / d
+ cy = r * cy / d
+ cz = r * cz / d
+ end if
+
+ end subroutine arc_bisect
+
+
+ subroutine poly_fit_2(a_in,b_out,weights_in,m,n,ne)
+
+ implicit none
+
+ integer, intent(in) :: m,n,ne
+ real (kind=RKIND), dimension(ne,ne), intent(in) :: a_in, weights_in
+ real (kind=RKIND), dimension(ne,ne), intent(out) :: b_out
+
+ ! local storage
+
+ real (kind=RKIND), dimension(m,n) :: a
+ real (kind=RKIND), dimension(n,m) :: b
+ real (kind=RKIND), dimension(m,m) :: w,wt,h
+ real (kind=RKIND), dimension(n,m) :: at, ath
+ real (kind=RKIND), dimension(n,n) :: ata, ata_inv, atha, atha_inv
+ integer, dimension(n) :: indx
+ integer :: i,j
+
+ if ( (ne<n) .or. (ne<m) ) then
+ write(6,*) ' error in poly_fit_2 inversion ',m,n,ne
+ stop
+ end if
+
+! a(1:m,1:n) = a_in(1:n,1:m)
+ a(1:m,1:n) = a_in(1:m,1:n)
+ w(1:m,1:m) = weights_in(1:m,1:m)
+ b_out(:,:) = 0.
+
+ wt = transpose(w)
+ h = matmul(wt,w)
+ at = transpose(a)
+ ath = matmul(at,h)
+ atha = matmul(ath,a)
+
+ ata = matmul(at,a)
+
+! if (m == n) then
+! call migs(a,n,b,indx)
+! else
+
+ call migs(atha,n,atha_inv,indx)
+
+ b = matmul(atha_inv,ath)
+
+! call migs(ata,n,ata_inv,indx)
+! b = matmul(ata_inv,at)
+! end if
+ b_out(1:n,1:m) = b(1:n,1:m)
+
+! do i=1,n
+! write(6,*) ' i, indx ',i,indx(i)
+! end do
+!
+! write(6,*) ' '
+
+ end subroutine poly_fit_2
+
+
+! Updated 10/24/2001.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! !
+! Please Note: !
+! !
+! (1) This computer program is written by Tao Pang in conjunction with !
+! his book, "An Introduction to Computational Physics," published !
+! by Cambridge University Press in 1997. !
+! !
+! (2) No warranties, express or implied, are made for this program. !
+! !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+SUBROUTINE MIGS (A,N,X,INDX)
+!
+! Subroutine to invert matrix A(N,N) with the inverse stored
+! in X(N,N) in the output. Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ INTEGER, INTENT (IN) :: N
+ INTEGER :: I,J,K
+ INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
+ REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
+ REAL (kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
+ REAL (kind=RKIND), DIMENSION (N,N) :: B
+!
+ DO I = 1, N
+ DO J = 1, N
+ B(I,J) = 0.0
+ END DO
+ END DO
+ DO I = 1, N
+ B(I,I) = 1.0
+ END DO
+!
+ CALL ELGS (A,N,INDX)
+!
+ DO I = 1, N-1
+ DO J = I+1, N
+ DO K = 1, N
+ B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K)
+ END DO
+ END DO
+ END DO
+!
+ DO I = 1, N
+ X(N,I) = B(INDX(N),I)/A(INDX(N),N)
+ DO J = N-1, 1, -1
+ X(J,I) = B(INDX(J),I)
+ DO K = J+1, N
+ X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
+ END DO
+ X(J,I) = X(J,I)/A(INDX(J),J)
+ END DO
+ END DO
+END SUBROUTINE MIGS
+
+
+SUBROUTINE ELGS (A,N,INDX)
+!
+! Subroutine to perform the partial-pivoting Gaussian elimination.
+! A(N,N) is the original matrix in the input and transformed matrix
+! plus the pivoting element ratios below the diagonal in the output.
+! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ INTEGER, INTENT (IN) :: N
+ INTEGER :: I,J,K,ITMP
+ INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
+ REAL (kind=RKIND) :: C1,PI,PI1,PJ
+ REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
+ REAL (kind=RKIND), DIMENSION (N) :: C
+!
+! Initialize the index
+!
+ DO I = 1, N
+ INDX(I) = I
+ END DO
+!
+! Find the rescaling factors, one from each row
+!
+ DO I = 1, N
+ C1= 0.0
+ DO J = 1, N
+ C1 = AMAX1(C1,ABS(A(I,J)))
+ END DO
+ C(I) = C1
+ END DO
+!
+! Search the pivoting (largest) element from each column
+!
+ DO J = 1, N-1
+ PI1 = 0.0
+ DO I = J, N
+ PI = ABS(A(INDX(I),J))/C(INDX(I))
+ IF (PI.GT.PI1) THEN
+ PI1 = PI
+ K = I
+ ENDIF
+ END DO
+!
+! Interchange the rows via INDX(N) to record pivoting order
+!
+ ITMP = INDX(J)
+ INDX(J) = INDX(K)
+ INDX(K) = ITMP
+ DO I = J+1, N
+ PJ = A(INDX(I),J)/A(INDX(J),J)
+!
+! Record pivoting ratios below the diagonal
+!
+ A(INDX(I),J) = PJ
+!
+! Modify other elements accordingly
+!
+ DO K = J+1, N
+ A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
+ END DO
+ END DO
+ END DO
+!
+END SUBROUTINE ELGS
+
+!-------------------------------------------------------------
+
+ subroutine initialize_deformation_weights( grid )
+
+!
+! compute the cell coefficients for the deformation calculations
+! WCS, 13 July 2010
+!
+ implicit none
+
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
+
+! local variables
+
+ real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
+ real (kind=RKIND), dimension(grid % nEdges) :: xe, ye
+ real (kind=RKIND), dimension(grid % nCells) :: theta_abs
+
+ real (kind=RKIND), dimension(25) :: xc, yc, zc ! cell center coordinates
+ real (kind=RKIND), dimension(25) :: thetav, thetat, dl_sphere
+ real (kind=RKIND) :: xm, ym, zm, dl, xec, yec, zec
+ real (kind=RKIND) :: thetae_tmp, xe_tmp, ye_tmp
+ real (kind=RKIND) :: xv1, xv2, yv1, yv2, zv1, zv2
+ integer :: i, j, k, ip1, ip2, m, n, ip1a, ii
+ integer :: iCell, iEdge
+ real (kind=RKIND) :: pii
+ real (kind=RKIND) :: x0, y0, x1, y1, x2, y2, x3, y3, x4, y4, x5, y5
+ real (kind=RKIND) :: pdx1, pdx2, pdx3, pdy1, pdy2, pdy3, dx1, dx2, dy1, dy2
+ real (kind=RKIND) :: angv1, angv2, dl1, dl2
+ real (kind=RKIND), dimension(25) :: dxe, dye, x2v, y2v, xp, yp, xpt, ypt
+
+ real (kind=RKIND) :: length_scale
+ integer :: ma,na, cell_add, mw, nn
+ integer, dimension(25) :: cell_list
+
+ integer :: cell1, cell2, iv
+ logical :: do_the_cell
+ real (kind=RKIND) :: area_cell, sint2, cost2, sint_cost, sumw1, sumw2, xptt, area_cellt
+
+ logical, parameter :: debug = .false.
+
+ if (debug) write(0,*) ' in def weight calc '
+
+ defc_a => grid % defc_a % array
+ defc_b => grid % defc_b % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+
+ defc_a(:,:) = 0.
+ defc_b(:,:) = 0.
+
+ pii = 2.*asin(1.0)
+
+ if (debug) write(0,*) ' beginning cell loop '
+
+ do iCell = 1, grid % nCells
+
+ if (debug) write(0,*) ' cell loop ', iCell
+
+ cell_list(1) = iCell
+ do i=2, grid % nEdgesOnCell % array(iCell)+1
+ cell_list(i) = grid % CellsOnCell % array(i-1,iCell)
+ end do
+ n = grid % nEdgesOnCell % array(iCell) + 1
+
+! check to see if we are reaching outside the halo
+
+ if (debug) write(0,*) ' points ', n
+
+ do_the_cell = .true.
+ do i=1,n
+ if (cell_list(i) > grid % nCells) do_the_cell = .false.
+ end do
+
+
+ if (.not. do_the_cell) cycle
+
+
+! 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
+
+
+ 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
+ end do
+
+ theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
+ xc(2), yc(2), zc(2), &
+ 0., 0., 1. )
+
+! angles from cell center to neighbor centers (thetav)
+
+ do i=1,n-1
+
+ ip2 = i+2
+ if (ip2 > n) ip2 = 2
+
+ thetav(i) = sphere_angle( xc(1), yc(1), zc(1), &
+ 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) )
+ end do
+
+ length_scale = 1.
+ do i=1,n-1
+ dl_sphere(i) = dl_sphere(i)/length_scale
+ end do
+
+ thetat(1) = 0. ! this defines the x direction, cell center 1 ->
+! thetat(1) = theta_abs(iCell) ! this defines the x direction, longitude line
+ do i=2,n-1
+ thetat(i) = thetat(i-1) + thetav(i-1)
+ end do
+
+ do i=1,n-1
+ xp(i) = cos(thetat(i)) * dl_sphere(i)
+ yp(i) = sin(thetat(i)) * dl_sphere(i)
+ end do
+
+ else ! On an x-y plane
+
+ xp(1) = grid % xCell % array(iCell)
+ yp(1) = grid % yCell % array(iCell)
+
+
+ do i=2,n
+ iv = grid % verticesOnCell % array(i-1,iCell)
+ xp(i) = grid % xVertex % array(iv)
+ yp(i) = grid % yVertex % array(iv)
+ end do
+
+ end if
+
+! thetat(1) = 0.
+ thetat(1) = theta_abs(iCell)
+ do i=2,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ thetat(i) = plane_angle( 0.,0.,0., &
+ xp(i)-xp(i-1), yp(i)-yp(i-1), 0., &
+ xp(ip1)-xp(i), yp(ip1)-yp(i), 0., &
+ 0., 0., 1.)
+ thetat(i) = thetat(i) + thetat(i-1)
+ end do
+
+ area_cell = 0.
+ area_cellt = 0.
+ do i=1,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
+ area_cell = area_cell + 0.25*(xp(i)+xp(ip1))*(yp(ip1)-yp(i)) - 0.25*(yp(i)+yp(ip1))*(xp(ip1)-xp(i))
+ area_cellt = area_cellt + (0.25*(xp(i)+xp(ip1))*cos(thetat(i)) + 0.25*(yp(i)+yp(ip1))*sin(thetat(i)))*dl
+ end do
+ if (debug) write(0,*) ' area_cell, area_cellt ',area_cell, area_cellt,area_cell-area_cellt
+
+ do i=1,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
+ sint2 = (sin(thetat(i)))**2
+ cost2 = (cos(thetat(i)))**2
+ sint_cost = sin(thetat(i))*cos(thetat(i))
+ defc_a(i,iCell) = dl*(cost2 - sint2)/area_cell
+ defc_b(i,iCell) = dl*2.*sint_cost/area_cell
+ if (cellsOnEdge(1,EdgesOnCell(i,iCell)) /= iCell) then
+ defc_a(i,iCell) = - defc_a(i,iCell)
+ defc_b(i,iCell) = - defc_b(i,iCell)
+ end if
+
+ end do
+
+ end do
+
+ if (debug) write(0,*) ' exiting def weight calc '
+
+ end subroutine initialize_deformation_weights
+
+end module advection
Added: branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_core.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_core.F (rev 0)
+++ branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_core.F 2010-11-03 18:27:30 UTC (rev 592)
@@ -0,0 +1,68 @@
+module core
+
+
+ contains
+
+
+ subroutine mpas_init(domain)
+
+ use grid_types
+ use configure
+ use test_cases
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ end subroutine mpas_init
+
+
+ subroutine mpas_run(domain, output_obj, output_frame)
+
+ use grid_types
+ use io_output
+ use timer
+ use test_cases
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ type (io_output_object), intent(inout) :: output_obj
+ integer, intent(inout) :: output_frame
+
+ integer :: ntimesteps, itimestep
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block
+
+
+ call setup_nhyd_test_case(domain)
+
+ !
+ ! Note: The following initialization calls have been moved to mpas_setup_test_case()
+ ! since values computed by these routines are needed to produce initial fields
+ !
+ ! call initialize_advection_rk(mesh)
+ ! call initialize_deformation_weights(mesh)
+
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % xtime % scalar = 0.0
+ block => block % next
+ end do
+
+ call output_state_for_domain(output_obj, domain, output_frame)
+
+ end subroutine mpas_run
+
+
+ subroutine mpas_finalize(domain)
+
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ end subroutine mpas_finalize
+
+end module core
Added: branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_test_cases.F (rev 0)
+++ branches/atmos_nonhydrostatic/src/core_init_nhyd_atmos/module_test_cases.F 2010-11-03 18:27:30 UTC (rev 592)
@@ -0,0 +1,1970 @@
+module test_cases
+
+ use grid_types
+ use configure
+ use constants
+ use dmpar
+ use advection
+
+
+ contains
+
+
+ subroutine setup_nhyd_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 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 nhyd_test_case_jw(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 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 nhyd_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 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 nhyd_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 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
+
+ end subroutine setup_nhyd_test_case
+
+!----------------------------------------------------------------------------------------------------------
+
+ subroutine nhyd_test_case_jw(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 :: 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 :: rh_max = 0.4 ! Maximum relative humidity
+ 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 :: zgrid, zx, zz, hx
+ real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt
+ real (kind=RKIND), dimension(:,:,:), pointer :: zf, zf3, zb, zb3
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+
+ 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, dimension(:), pointer :: dvEdge, AreaCell
+ real, 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), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature, qv
+ real (kind=RKIND) :: ptmp, 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
+
+ real (kind=RKIND) :: d1, d2, d3, cf1, cf2, cf3, cof1, cof2, psurf
+
+ ! storage for (lat,z) arrays for zonal velocity calculation
+
+ integer, parameter :: nlat=361
+ real (kind=RKIND), dimension(grid % nVertLevels + 1) :: zz_1d, zgrid_1d, hx_1d
+ real (kind=RKIND), dimension(grid % nVertLevels) :: flux_zonal
+ real (kind=RKIND), dimension(nlat, grid % nVertLevels) :: u_2d, etavs_2d
+ real (kind=RKIND), dimension(nlat) :: lat_2d
+ real (kind=RKIND) :: dlat
+ real (kind=RKIND) :: z_edge, z_edge3, d2fdx2_cell1, d2fdx2_cell2
+
+ !
+ ! 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
+ zf => grid % zf % array
+ zf3 => grid % zf3% 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 => state % rho % array
+ rr => diag % rho_p % array
+ t => state % theta % array
+ rt => diag % rtheta_p % array
+
+ scalars => state % scalars % array
+
+ scalars(:,:,:) = 0.
+
+ call initialize_advection_rk(grid)
+ call initialize_deformation_weights(grid)
+
+ 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.5
+ 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 (lat,z) calc for zonal velocity
+
+ dlat = 0.5*pii/float(nlat-1)
+ do i = 1,nlat
+
+ lat_2d(i) = float(i-1)*dlat
+! write(0,*) ' zonal setup, latitude = ',lat_2d(i)*180./pii
+
+ do k=1,nz
+ phi = lat_2d(i)
+ hx_1d(k) = 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
+
+ do k=1,nz
+ zgrid_1d(k) = (1.-ah(k))*(sh(k)*(zt-hx_1d(k))+hx_1d(k)) &
+ + ah(k) * sh(k)* zt
+ end do
+ do k=1,nz1
+ zz_1d (k) = (zw(k+1)-zw(k))/(zgrid_1d(k+1)-zgrid_1d(k))
+ end do
+
+ do k=1,nz1
+ ztemp = .5*(zgrid_1d(k+1)+zgrid_1d(k))
+ 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_1d(k))
+ 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 = grid % latCell % array (i)
+ phi = lat_2d (i)
+ do k=1,nz1
+ tt(k) = 0.
+ tt(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)
+
+
+ ztemp = .5*(zgrid_1d(k)+zgrid_1d(k+1))
+ ptemp = ppb(k,i) + pp(k,i)
+ qv(k,i) = 0.
+
+ end do
+
+ do itrp = 1,25
+ do k=1,nz1
+ rr(k,i) = (pp(k,i)/(rgas*zz_1d(k)) &
+ -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(1,i)) &
+ -.25*(rr(2,i)+rb(2,i))*(1.+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))*qv(k ,i) &
+ +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*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
+ etavs_2d(i,k) = (0.5*(ppb(k,i)+ppb(k,i)+pp(k,i)+pp(k,i))/p0 - 0.252)*pii/2.
+! u_2d(i,k) = u0*(sin(2.*lat_2d(i))**2) *(cos(etavs_2d(i,k))**1.5)
+ u_2d(i,k) = u0*(sin(2.*lat_2d(i))**2) *(cos(etavs_2d(i,k))**1.5)*(rb(k,i)+rr(k,i))
+ end do
+
+ end do ! end loop over latitudes for 2D zonal wind field calc
+
+! do i=1,nlat
+! do k=1,nz1
+! u_2d(i,k) = u_2d(i,k) - u0*(sin(2.*lat_2d(i))**2) *(cos(etavs_2d(nlat/2,k))**1.5)
+! end do
+! end do
+!
+! write(22,*) nz1,nlat,u_2d
+
+!******************************************************************
+
+!
+!---- baroclinc wave initialization ---------------------------------
+!
+! reference sounding based on dry isothermal atmosphere
+!
+ do i=1, grid % nCells
+ !write(0,*) ' thermodynamic setup, cell ',i
+ 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
+!
+! 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
+ tt(k) = 0.
+ tt(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)
+
+
+ !write(0,*) ' k, tt(k) ',k,tt(k)
+ ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
+ ptemp = ppb(k,i) + pp(k,i)
+! qv(k,i) = env_qv( ztemp, tt(k), ptemp, 0 )
+ qv(k,i) = 0.
+
+ end do
+! do k=2,nz1
+! cqw(k,i) = 1./(1.+.5*(qv(k,i)+qv(k-1,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.+qv(1,i)) &
+ -.25*(rr(2,i)+rb(2,i))*(1.+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))*qv(k ,i) &
+ +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*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 (k,i) = rb(k,i) + rr(k,i)
+ end do
+
+ if(i == 1) then
+ do k=1,nz1
+ write(0,*) ' k, p, t, rt ',k,p(k,1),t(k,1),rt(k,1)
+ enddo
+ end if
+
+ end do ! end loop over cells
+
+ 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.)/(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
+
+ call 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
+!! etavs = (0.5*(ppb(k,iCell1)+ppb(k,iCell2)+pp(k,iCell1)+pp(k,iCell2))/p0 - 0.252)*pii/2.
+! etavs = (0.5*(ppb(k,1)+ppb(k,1)+pp(k,1)+pp(k,1))/p0 - 0.252)*pii/2.
+ etavs = (0.5*(ppb(k,440)+ppb(k,440)+pp(k,440)+pp(k,440))/p0 - 0.252)*pii/2. ! 10262 mesh
+! etavs = (0.5*(ppb(k,505)+ppb(k,505)+pp(k,505)+pp(k,505))/p0 - 0.252)*pii/2. ! 40962 mesh
+
+! fluxk = u0*flux*(cos(etavs)**1.5)
+
+ fluxk = u0*flux_zonal(k)/(0.5*(rb(k,iCell1)+rb(k,iCell2)+rr(k,iCell1)+rr(k,iCell2)))
+
+! if(k.eq.18) then
+! write(21,*) ' iEdge, u1, u2 ',iEdge,fluxk,u0*flux_zonal(k)
+! end if
+!! fluxk = u0*flux*(cos(znuv(k))**(1.5))
+!! fluxk = u0 * cos(grid % angleEdge % array(iEdge)) * (sin(lat1+lat2)**2) *(cos(etavs)**1.5)
+ state % u % array(k,iEdge) = fluxk + u_pert
+ end do
+
+ cell1 = grid % CellsOnEdge % array(1,i)
+ cell2 = grid % CellsOnEdge % array(2,i)
+ if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+ do k=1,nz1
+ diag % ru % array (k,iEdge) = 0.5*(rho(k,cell1)+rho(k,cell2))*state % u % array (k,iEdge)
+ end do
+ end if
+
+ !
+ ! 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)
+ 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 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)
+ 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
+ diag % rw % array = 0.
+ state % w % array = 0.
+ 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)*diag % ru % array(k,iEdge)+fzp(k)*diag % ru % array(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.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
+ + sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ end if
+
+ end do
+ end if
+
+ end do
+
+ ! Compute w from rho 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 % array(k-1,iCell) + fzm(k) * state % rho % 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)
+ 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 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.+qv(1,i)) &
+ -.25*(rr(2,i)+rb(2,i))*(1.+qv(2,i)))
+
+ write(0,*) ' i, psurf, lat ',i,psurf,grid%latCell%array(i)*180./3.1415828
+ enddo
+
+ end subroutine nhyd_test_case_jw
+
+ subroutine 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(nlat,nz1), 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(i,k) + w2*u_2d(i+1,k)
+ end do
+
+ end if
+
+ end do
+
+! renormalize for setting cell-face fluxes
+
+ do k=1,nz1
+ flux_zonal(k) = sign(1.,lat2_in-lat1_in)*flux_zonal(k)*dlat*a/dvEdge/u0
+ end do
+
+ end subroutine calc_flux_zonal
+
+
+!----------------------------------------------------------------------------------------------------------
+
+ subroutine nhyd_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, 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, 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 => state % rho % array
+
+ pp => diag % pressure_p % array
+ rr => diag % rho_p % array
+ t => state % theta % 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 initialize_advection_rk(grid)
+ call 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 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 dmpar_bcast_real(dminfo, pitop)
+ call 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)
+ 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) = amin1(0.014,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 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 dmpar_max_real(dminfo, maxval(grid % xCell % array(:)), xmid)
+ call 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 dmpar_bcast_real(dminfo, pitop)
+
+ 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
+ print *, "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 dmpar_bcast_reals(dminfo, nz1, t_init_1d)
+ call 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(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(k,cell1)+rho(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 % zf % array = 0.
+ grid % zf3% 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
+
+ end subroutine nhyd_test_case_squall_line
+
+
+!----------------------------------------------------------------------------------------------------------
+
+
+ subroutine nhyd_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, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
+ real (kind=RKIND), dimension(:,:,:), pointer :: scalars, deriv_two, zf, zf3, 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, dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell
+ real, 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) :: ptmp, 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
+ zf => grid % zf % array
+ zf3 => grid % zf3 % 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 => state % rho % array
+
+ pp => diag % pressure_p % array
+ rr => diag % rho_p % array
+ t => state % theta % 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 initialize_advection_rk(grid)
+ call 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.)
+ 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 print *,"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
+ u(k,i) = cos(grid % angleEdge % array(i)) * (u(k,i) - us)
+ 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) = amin1(0.014,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(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(k,cell1)+rho(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)
+
+ 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.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
+ + sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ end if
+
+ end do
+ end if
+
+ end do
+
+ ! Compute w from rho 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 % array(k-1,iCell) + fzm(k) * state % rho % 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
+
+ end subroutine nhyd_test_case_mtn_wave
+
+
+!----------------------------------------------------------------------------------------------------------
+
+ real 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
+
+end module test_cases
</font>
</pre>