<p><b>laura@ucar.edu</b> 2010-08-03 12:39:56 -0600 (Tue, 03 Aug 2010)</p><p>Merged with nonhydrostatic branch - revision 444<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/Registry 2010-08-03 18:36:36 UTC (rev 457)
+++ branches/atmos_physics/src/core_nhyd_atmos/Registry 2010-08-03 18:39:56 UTC (rev 458)
@@ -1,25 +1,35 @@
#
# namelist type namelist_record name default_value
#
-namelist integer sw_model config_test_case 5
-namelist character sw_model config_time_integration SRK3
-namelist real sw_model config_dt 172.8
-namelist integer sw_model config_ntimesteps 7500
-namelist integer sw_model config_output_interval 500
-namelist real sw_model config_h_mom_eddy_visc2 0.0
-namelist real sw_model config_h_mom_eddy_visc4 0.0
-namelist real sw_model config_v_mom_eddy_visc2 0.0
-namelist real sw_model config_h_theta_eddy_visc2 0.0
-namelist real sw_model config_h_theta_eddy_visc4 0.0
-namelist real sw_model config_v_theta_eddy_visc2 0.0
-namelist integer sw_model config_number_of_sub_steps 4
-namelist integer sw_model config_theta_adv_order 2
-namelist integer sw_model config_scalar_adv_order 2
-namelist logical sw_model config_positive_definite false
-namelist logical sw_model config_monotonic true
-namelist integer sw_model config_mp_physics 0
-namelist real sw_model config_epssm 0.1
-namelist real sw_model config_smdiv 0.1
+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 character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -225,6 +235,11 @@
var real deriv_two ( FIFTEEN TWO nEdges ) o deriv_two - -
var integer advCells ( TWENTYONE nCells ) - advCells - -
+# Space needed for deformation calculation weights
+var real defc_a ( maxEdges nCells ) - defc_a - -
+var real defc_b ( maxEdges nCells ) - defc_b - -
+var real kdiff ( nVertLevels nCells Time ) - kdiff - -
+
# Arrays required for reconstruction of velocity field
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2010-08-03 18:36:36 UTC (rev 457)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2010-08-03 18:39:56 UTC (rev 458)
@@ -112,7 +112,7 @@
real (kind=RKIND) :: r_earth, etavs, ztemp, zd, zt, dz, gam, delt, str
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh, temperature, qv
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature, qv
real (kind=RKIND) :: ptmp, es, qvs, xnutr, znut, ptemp
integer :: iter
@@ -134,8 +134,6 @@
real (kind=RKIND), dimension(nlat) :: lat_2d
real (kind=RKIND) :: dlat
- integer:: iq
-
!
! Scale all distances and areas from a unit sphere to one with radius a
!
@@ -187,6 +185,7 @@
rt => grid % rtheta_p % array
scalars => state % scalars % array
+
scalars(:,:,:) = 0.
xnutr = 0.
@@ -360,10 +359,6 @@
end do
-!
-! iterations to converge temperature as a function of pressure
-!
-
do itr = 1,10
do k=1,nz1
@@ -390,6 +385,7 @@
ztemp = .5*(zgrid_1d(k)+zgrid_1d(k+1))
ptemp = ppb(k,i) + pp(k,i)
+ qv(k,i) = 0.
end do
@@ -435,8 +431,7 @@
! write(22,*) nz1,nlat,u_2d
!******************************************************************
- write(6,*)
- write(6,*) '--- RELATIVE HUMIDITY:', zt
+
!
!---- baroclinc wave initialization ---------------------------------
!
@@ -456,11 +451,11 @@
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
+ 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
#ifdef DO_PHYSICS
!
@@ -470,44 +465,21 @@
201 format(3i6,8(1x,e15.8))
202 format(2i6,8(1x,e15.8))
- if(config_microp_scheme .ne. 'off' .or. &
- config_conv_shallow_scheme .ne. 'off' .or. &
- config_conv_deep_scheme .ne. 'off' .or. &
- config_pbl_scheme .ne. 'off' .or. &
- config_eddy_scheme .ne. 'off' .or. &
- config_radt_lw_scheme .ne. 'off' .or. &
- config_radt_sw_scheme .ne. 'off') then
-
-!... relative humidity as function of pressure:
-! do k = 1, nz1
-! ptemp = ppb(k,i) + pp(k,i)
-! if(ptemp < 50000.) then
-! rh(k,i) = 0.0
-! elseif(p0-ptemp .lt. 0.) then
-! rh(k,i) = 1.0
-! else
-! rh(k,i) = (1.-((p0-ptemp)/50000.)**1.25)
-! endif
-! rh(k,i) = min(rh_max,rh(k,i))
-! enddo
-
-!... relative humidity as function of height:
- do k = 1, nz1
- ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
- if(ztemp .gt. 5000.) then
- rh(k,i) = 0.0
- else
- rh(k,i) = (1.-0.75*(ztemp/zt)**1.25)
- endif
- rh(k,i) = min(rh_max,rh(k,i))
-! write(6,202) i,k,ztemp,rh(k,i)
- enddo
- endif
+ do k = 1, nz1
+ ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
+ if(ztemp .gt. 5000.) then
+ rel_hum(k,i) = 0.0
+ elseif(ztemp .lt. 0.0) then
+ rel_hum(k,i) = 1.0
+ else
+ rel_hum(k,i) = (1.-0.75*(ztemp/zt)**1.25)
+ endif
+ rel_hum(k,i) = min(rh_max,rel_hum(k,i))
+ enddo
#endif
!
! iterations to converge temperature as a function of pressure
!
-! write(6,*)
do itr = 1,10
do k=1,nz1
@@ -534,6 +506,8 @@
!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.
#ifdef DO_PHYSICS
!
@@ -546,11 +520,10 @@
es = 1000.*0.6112*exp(21.8745584*(tt(k)-273.16)/(tt(k)-7.66))
end if
qvs = (287.04/461.6)*es/(ptemp-es)
- scalars(index_qv,k,i) = rh(k,i)*qvs
-! write(6,201) itr,i,k,ztemp,ptemp,es,qvs,rh(k,i),scalars(index_qv,k,i)
+ scalars(index_qv,k,i) = rel_hum(k,i)*qvs
+! write(6,201) itr,i,k,ztemp,zgrid(k,i),zgrid(k+1,i),ptemp,tt(k),qvs,rel_hum(k,i),scalars(index_qv,k,i)
#endif
end do
-! write(0,*)
! do k=2,nz1
! cqw(k,i) = 1./(1.+.5*(qv(k,i)+qv(k-1,i)))
@@ -588,11 +561,11 @@
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
+ 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
@@ -660,6 +633,15 @@
end do
!
+ ! CALCULATION OF OMEGA, RW = ZX * RU + ZZ * RW
+ !
+ ! we are assuming w and rw are zero for this initialization
+ ! i.e., no terrain
+ !
+ grid % rw % array = 0.
+ state % w % array = 0.
+
+ !
! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
!
state % v % array(:,:) = 0.0
@@ -685,14 +667,6 @@
enddo
! stop
-! write(6,*)
-! write(6,*) '--- END JW:'
-! do iCell = 1, grid%nCells
-! do k = 1,grid % nVertLevels
-! write(6,202) iCell,k,(state%scalars%array(iq,k,iCell),iq=moist_start,moist_end)
-! enddo
-! 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)
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2010-08-03 18:36:36 UTC (rev 457)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2010-08-03 18:39:56 UTC (rev 458)
@@ -79,9 +79,8 @@
logical, parameter :: debug = .false.
! logical, parameter :: debug = .true.
logical, parameter :: debug_mass_conservation = .true.
- logical, parameter :: do_microphysics = .false.
- logical, parameter :: scalar_advection = .false.
-
+ logical, parameter :: do_microphysics = .false.
+ logical, parameter :: scalar_advection = .true.
real (kind=RKIND) :: domain_mass, scalar_mass, scalar_min, scalar_max
real (kind=RKIND) :: global_domain_mass, global_scalar_mass, global_scalar_min, global_scalar_max
@@ -254,7 +253,7 @@
! ************ advection of moist variables here...
- if(scalar_advection) then
+ if (config_scalar_advection) then
block => domain % blocklist
do while (associated(block))
@@ -794,7 +793,7 @@
- rdzw(k)*(dpzx(k+1)-dpzx(k))
pgrad = 0.5*c2*(exner(k,cell1)+exner(k,cell2))*pgrad
du(k) = dts*(tend_ru(k,iEdge) - cqu(k,iEdge) * pgrad)
-! + (0.005/6.)*dcEdge(iEdge)*(h_divergence(k,cell2)-h_divergence(k,cell1))
+! + (0.05/6.)*dcEdge(iEdge)*(h_divergence(k,cell2)-h_divergence(k,cell1))
ru_p(k,iEdge) = ru_p(k,iEdge) + du(k)
@@ -1094,7 +1093,7 @@
real (kind=RKIND), dimension(:,:,:), pointer :: scalar_old, scalar_new, scalar_tend
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND), dimension(:,:), pointer :: uhAvg, h_old, h_new, wwAvg, rho_edge, rho, zgrid
+ real (kind=RKIND), dimension(:,:), pointer :: uhAvg, h_old, h_new, wwAvg, rho_edge, rho, zgrid, kdiff
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, qv_init
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -1104,17 +1103,27 @@
real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw
real (kind=RKIND) :: coef_3rd_order
-
real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2, scalar_turb_flux, z1,z2,z3,z4,zm,z0,zp
- logical, parameter :: mix_full = .false.
-! logical, parameter :: mix_full = .true.
- coef_3rd_order = 0.
- if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
- if (config_scalar_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+ real (kind=RKIND) :: flux3, flux4
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
+ ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+
+ flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
+ coef3*ua*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+
+! coef_3rd_order = 0.
+! if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
+! if (config_scalar_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
+ coef_3rd_order = config_coef_3rd_order
+
scalar_old => s_old % scalars % array
scalar_new => s_new % scalars % array
+ kdiff => s_new % kdiff % array
deriv_two => grid % deriv_two % array
!**** uhAvg => grid % uhAvg % array
uhAvg => grid % ruAvg % array
@@ -1205,17 +1214,6 @@
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
end if
-! old version of the above code, with coef_3rd_order assumed to be 1.0
-! if (uhAvg(k,iEdge) > 0) then
-! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
-! 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-! -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
-! else
-! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
-! 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-! -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
-! end if
-
scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) - flux/areaCell(cell1)
scalar_tend(iScalar,k,cell2) = scalar_tend(iScalar,k,cell2) + flux/areaCell(cell2)
@@ -1282,6 +1280,26 @@
end if
end do
+ else if ( config_horiz_mixing == "2d_smagorinsky") then
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,grid % nVertLevels
+ do iScalar=1,num_scalars
+ scalar_turb_flux = 0.5*(kdiff(k,cell1)+kdiff(k,cell2))*prandtl* &
+ (scalar_new(iScalar,k,cell2) - scalar_new(iScalar,k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * rho_edge(k,iEdge) * scalar_turb_flux
+ scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) + flux/areaCell(cell1)
+ scalar_tend(iScalar,k,cell2) = scalar_tend(iScalar,k,cell2) - flux/areaCell(cell2)
+ end do
+ end do
+
+ end if
+ end do
+
end if
! vertical mixing
@@ -1306,7 +1324,7 @@
end do
end do
- if ( .not. mix_full) then
+ if ( .not. config_mix_full) then
iScalar = index_qv
do k=2,nVertLevels-1
z1 = zgrid(k-1,iCell)
@@ -1334,14 +1352,42 @@
do iCell=1,grid % nCells
- wdtn(:,1) = 0.
- do k = 2, nVertLevels
- do iScalar=1,num_scalars
- wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
- end do
- end do
- wdtn(:,nVertLevels+1) = 0.
+ wdtn(:,1) = 0.
+ if (config_scalar_vadv_order == 2) then
+ do k = 2, nVertLevels
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ end do
+ end do
+ else
+ k = 2
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ enddo
+ if ( config_scalar_vadv_order == 3 ) then
+ do k=3,nVertLevels-1
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = flux3( scalar_new(iScalar,k-2,iCell),scalar_new(iScalar,k-1,iCell), &
+ scalar_new(iScalar,k ,iCell),scalar_new(iScalar,k+1,iCell), &
+ wwAvg(k,iCell), coef_3rd_order )
+ end do
+ end do
+ else if ( config_scalar_vadv_order == 4 ) then
+ do k=3,nVertLevels-1
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = flux4( scalar_new(iScalar,k-2,iCell),scalar_new(iScalar,k-1,iCell), &
+ scalar_new(iScalar,k ,iCell),scalar_new(iScalar,k+1,iCell), wwAvg(k,iCell) )
+ end do
+ end do
+ end if
+ k = nVertLevels
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ enddo
+ end if
+ wdtn(:,nVertLevels+1) = 0.
+
do k=1,grid % nVertLevelsSolve
do iScalar=1,num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
@@ -1395,6 +1441,18 @@
real (kind=RKIND), parameter :: eps=1.e-20
real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: flux3, flux4
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
+ ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+
+ flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
+ coef3*ua*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+
+!------
+
scalar_old => s_old % scalars % array
scalar_new => s_new % scalars % array
deriv_two => grid % deriv_two % array
@@ -1433,10 +1491,12 @@
scale_out(:,:,:) = 1.
scale_in(:,:,:) = 1.
- coef_3rd_order = 0.
- if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
- if (config_scalar_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+! coef_3rd_order = 0.
+! if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
+! if (config_scalar_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+ coef_3rd_order = config_coef_3rd_order
+
do k = 1, grid % nVertLevels
kcp1 = min(k+1,grid % nVertLevels)
kcm1 = max(k-1,1)
@@ -1446,17 +1506,37 @@
do iCell=1,grid % nCells
if (k < grid % nVertLevels) then
+
+ if ((config_scalar_vadv_order == 2) .or. (k==1) .or. (k==grid % nVertLevels)) then
+ do iScalar=1,num_scalars
+ v_flux(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) * &
+ (fnm(k+1) * scalar_new(iScalar,k+1,iCell) + fnp(k+1) * scalar_new(iScalar,k,iCell))
+ end do
+ else if (config_scalar_vadv_order == 3) then
+ do iScalar=1,num_scalars
+ v_flux(iScalar,iCell,km0) = dt * flux3( scalar_new(iScalar,k-1,iCell),scalar_new(iScalar,k ,iCell), &
+ scalar_new(iScalar,k+1,iCell),scalar_new(iScalar,k+2,iCell), &
+ wwAvg(k+1,iCell), coef_3rd_order )
+ end do
+ else if (config_scalar_vadv_order == 4) then
+ do iScalar=1,num_scalars
+ v_flux(iScalar,iCell,km0) = dt * flux4( scalar_new(iScalar,k-1,iCell),scalar_new(iScalar,k ,iCell), &
+ scalar_new(iScalar,k+1,iCell),scalar_new(iScalar,k+2,iCell), wwAvg(k+1,iCell) )
+ end do
+ end if
+
cell_upwind = k
if (wwAvg(k+1,iCell) >= 0) cell_upwind = k+1
+
do iScalar=1,num_scalars
- v_flux(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) * &
- (fnm(k+1) * scalar_new(iScalar,k+1,iCell) + fnp(k+1) * scalar_new(iScalar,k,iCell))
v_flux_upwind(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) * scalar_old(iScalar,cell_upwind,iCell)
v_flux(iScalar,iCell,km0) = v_flux(iScalar,iCell,km0) - v_flux_upwind(iScalar,iCell,km0)
! v_flux(iScalar,iCell,km0) = 0. ! use only upwind - for testing
s_update(iScalar,iCell,km0) = scalar_old(iScalar,k,iCell) * h_old(k,iCell) &
- rdnw(k) * (v_flux_upwind(iScalar,iCell,km0) - v_flux_upwind(iScalar,iCell,km1))
end do
+
+
else
do iScalar=1,num_scalars
v_flux(iScalar,iCell,km0) = 0.
@@ -1733,7 +1813,7 @@
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, zgrid, rho_edge, rho, ru, u, v, tend_u, &
circulation, divergence, vorticity, ke, pv_edge, theta, rw, tend_rho, &
rt_diabatic, tend_theta, tend_w, w, cqw, rb, rr, pp, pressure_b, zz, zx, cqu, &
- h_divergence
+ h_divergence, kdiff
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -1753,17 +1833,36 @@
! logical, parameter :: debug = .true.
logical, parameter :: debug = .false.
- logical, parameter :: mix_full = .false.
-! logical, parameter :: mix_full = .true.
- integer :: w_adv_order
+ real (kind=RKIND), parameter :: c_s = 0.25
+ real (kind=RKIND), dimension( grid % nVertLevels ) :: d_diag, d_off_diag
+ real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
+ logical :: delsq_horiz_mixing
+
+
real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: flux3, flux4
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
+ ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+
+ flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
+ coef3*ua*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+
+!-----------
+
+ coef_3rd_order = config_coef_3rd_order
+
rho => s % rho % array
rho_edge => s % rho_edge % array
rb => grid % rho_base % array
rr => s % rho_p % array
u => s % u % array
+ v => s % v % array
+ kdiff => s % kdiff % array
ru => grid % ru % array
w => s % w % array
rw => grid % rw % array
@@ -1783,6 +1882,7 @@
verticesOnEdge => grid % verticesOnEdge % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
dcEdge => grid % dcEdge % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
@@ -1792,6 +1892,10 @@
zz => grid % zz % array
zx => grid % zx % array
+ defc_a => grid % defc_a % array
+ defc_b => grid % defc_b % array
+
+
tend_u => tend % u % array
tend_theta => tend % theta % array
tend_w => tend % w % array
@@ -1863,6 +1967,28 @@
end do
end do
+ delsq_horiz_mixing = .false.
+ if (config_horiz_mixing == "2d_smagorinsky") then
+ do iCell = 1, nCells
+ d_diag(:) = 0.
+ d_off_diag(:) = 0.
+ do iEdge = 1, grid % nEdgesOnCell % array (iCell)
+ do k=1, nVertLevels
+ d_diag(k) = d_diag(k) + defc_a(iEdge,iCell)*u(k,EdgesOnCell(iEdge,iCell)) &
+ - defc_b(iEdge,iCell)*v(k,EdgesOnCell(iEdge,iCell))
+ d_off_diag(k) = d_off_diag(k) + defc_b(iEdge,iCell)*u(k,EdgesOnCell(iEdge,iCell)) &
+ + defc_a(iEdge,iCell)*v(k,EdgesOnCell(iEdge,iCell))
+ end do
+ end do
+ do k=1, nVertLevels
+ kdiff(k,iCell) = (c_s * config_len_disp)**2 * sqrt(d_diag(k)**2 + d_off_diag(k)**2)
+ end do
+ end do
+ delsq_horiz_mixing = .true.
+ else if ( config_horiz_mixing == "2d_fixed") then
+ delsq_horiz_mixing = .true.
+ end if
+
#ifdef LANL_FORMULATION
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
@@ -1959,6 +2085,37 @@
cell2 = cellsOnEdge(2,iEdge)
wduz(1) = 0.
+ if (config_u_vadv_order == 2) then
+
+ do k=2,nVertLevels
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+ end do
+
+ else if (config_u_vadv_order == 3) then
+
+ k = 2
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+ do k=3,nVertLevels-1
+ wduz(k) = flux3( u(k-2,iEdge),u(k-1,iEdge),u(k,iEdge),u(k+1,iEdge),0.5*(rw(k,cell1)+rw(k,cell2)), 1. )
+ end do
+ k = nVertLevels
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+
+ else if (config_u_vadv_order == 4) then
+
+ k = 2
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+ do k=3,nVertLevels-1
+ wduz(k) = flux4( u(k-2,iEdge),u(k-1,iEdge),u(k,iEdge),u(k+1,iEdge),0.5*(rw(k,cell1)+rw(k,cell2)) )
+ end do
+ k = nVertLevels
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+
+ end if
+
+ wduz(nVertLevels+1) = 0.
+
+ wduz(1) = 0.
do k=2,nVertLevels
wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
end do
@@ -1972,28 +2129,54 @@
!
! horizontal mixing for u
!
- if ( h_mom_eddy_visc2 > 0.0 ) then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ if (delsq_horiz_mixing) then
- !
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
- ! only valid for h_mom_eddy_visc2 == constant
- !
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = rho_edge(k,iEdge)*h_mom_eddy_visc2 * u_diffusion
+ if (h_mom_eddy_visc2 > 0.0) then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
+ ! only valid for h_mom_eddy_visc2 == constant
+ !
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = rho_edge(k,iEdge)*h_mom_eddy_visc2 * u_diffusion
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ end do
+ end do
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
- end do
- end do
- end if
+ else if ( config_horiz_mixing == "2d_smagorinsky") then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
+ ! only valid for h_mom_eddy_visc2 == constant
+ !
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = rho_edge(k,iEdge)* 0.5*(kdiff(k,cell1)+kdiff(k,cell2)) * u_diffusion
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ end do
+ end do
+ end if
+
+ end if ! delsq_horiz_mixing for u
+
if ( h_mom_eddy_visc4 > 0.0 ) then
allocate(delsq_divergence(nVertLevels, nCells+1))
@@ -2097,7 +2280,7 @@
!
if ( v_mom_eddy_visc2 > 0.0 ) then
- if (mix_full) then
+ if (config_mix_full) then
do iEdge=1,grid % nEdgesSolve
@@ -2161,10 +2344,8 @@
! horizontal advection for w
!
- w_adv_order = 2
+ if (config_w_adv_order == 2) then
- if (w_adv_order == 2) then
-
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -2178,7 +2359,7 @@
end if
end do
- else if (w_adv_order == 3) then
+ else if (config_w_adv_order == 3) then
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -2216,7 +2397,7 @@
end if
end do
- else if (w_adv_order == 4) then
+ else if (config_w_adv_order == 4) then
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -2256,9 +2437,11 @@
! Note: we are using quite a bit of the theta code here - could be combined later???
- if ( h_mom_eddy_visc2 > 0.0 ) then
+ if (delsq_horiz_mixing) then
- do iEdge=1,grid % nEdges
+ if (h_mom_eddy_visc2 > 0.0) then
+
+ do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
@@ -2270,8 +2453,27 @@
tend_w(k,cell2) = tend_w(k,cell2) - flux
end do
- end if
- end do
+ end if
+ end do
+ else if (config_horiz_mixing == "2d_smagorinsky") then
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+
+ do k=2,grid % nVertLevels
+! theta_turb_flux = h_mom_eddy_visc2*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = 0.25*(kdiff(k,cell1)+kdiff(k,cell2)+kdiff(k-1,cell1)+kdiff(k,cell2)) &
+ *(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
+ flux = 0.5*dvEdge (iEdge) * (rho_edge(k,iEdge)+rho_edge(k-1,iEdge)) * theta_turb_flux
+ tend_w(k,cell1) = tend_w(k,cell1) + flux
+ tend_w(k,cell2) = tend_w(k,cell2) - flux
+ end do
+
+ end if
+ end do
+ end if
end if
@@ -2327,14 +2529,40 @@
!
do iCell = 1, nCells
+
wdwz(1) = 0.
- do k=2,nVertLevels
+ if (config_w_vadv_order == 2) then
+
+ do k=2,nVertLevels
+ wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
+ end do
+
+ else if (config_w_vadv_order == 3) then
+
+ k = 2
wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
- end do
+ do k=3,nVertLevels-1
+ wdwz(k) = flux3( w(k-2,iCell),w(k-1,iCell),w(k,iCell),w(k+1,iCell),0.5*(rw(k,iCell)+rw(k-1,iCell)), 1. )
+ end do
+ k = nVertLevels
+ wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
+
+ else if (config_w_vadv_order == 4) then
+
+ k = 2
+ wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
+ do k=3,nVertLevels-1
+ wdwz(k) = flux4( w(k-2,iCell),w(k-1,iCell),w(k,iCell),w(k+1,iCell),0.5*(rw(k,iCell)+rw(k-1,iCell)) )
+ end do
+ k = nVertLevels
+ wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
+
+ end if
+
wdwz(nVertLevels+1) = 0.
+
do k=2,nVertLevels
-
tend_w(k,iCell) = tend_w(k,iCell)/areaCell(iCell) -rdzu(k)*(wdwz(k+1)-wdwz(k)) &
- cqw(k,iCell)*( rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
+ gravity* &
@@ -2424,8 +2652,6 @@
! 3rd order stencil
- coef_3rd_order = 0.25
-
if( u(k,iEdge) > 0) then
flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
0.5*(theta(k,cell1) + theta(k,cell2)) &
@@ -2491,22 +2717,42 @@
! horizontal mixing for theta - we could combine this with advection directly (i.e. as a turbulent flux),
! but here we can also code in hyperdiffusion if we wish (2nd order at present)
!
- if ( h_theta_eddy_visc2 > 0.0 ) then
+ if (delsq_horiz_mixing) then
+ if ( h_theta_eddy_visc2 > 0.0 ) then
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+
+ do k=1,grid % nVertLevels
+ theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
+ tend_theta(k,cell1) = tend_theta(k,cell1) + flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+ end do
+
+ end if
+ end do
- do k=1,grid % nVertLevels
- theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
- flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
- tend_theta(k,cell1) = tend_theta(k,cell1) + flux
- tend_theta(k,cell2) = tend_theta(k,cell2) - flux
- end do
+ else if ( ( config_horiz_mixing == "2d_smagorinsky") ) then
- end if
- end do
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+
+ do k=1,grid % nVertLevels
+ theta_turb_flux = 0.5*(kdiff(k,cell1)+kdiff(k,cell2))*prandtl &
+ *(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
+ tend_theta(k,cell1) = tend_theta(k,cell1) + flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+ end do
+
+ end if
+ end do
+ end if
end if
@@ -2561,11 +2807,39 @@
! Note: we are also dividing through by the cell area after the horizontal flux divergence
!
do iCell = 1, nCells
+
wdtz(1) = 0.
- do k=2,nVertLevels
+ if (config_theta_vadv_order == 2) then
+
+ do k=2,nVertLevels
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
+ end do
+
+ else if (config_theta_vadv_order == 3) then
+
+ k = 2
wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
- end do
+ do k=3,nVertLevels-1
+ wdtz(k) = flux3( theta(k-2,iCell),theta(k-1,iCell),theta(k,iCell),theta(k+1,iCell), rw(k,iCell), coef_3rd_order )
+ end do
+ k = nVertLevels
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
+
+ else if (config_theta_vadv_order == 4) then
+
+ k = 2
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
+ do k=3,nVertLevels-1
+ wdtz(k) = flux4( theta(k-2,iCell),theta(k-1,iCell),theta(k,iCell),theta(k+1,iCell), rw(k,iCell) )
+ end do
+ k = nVertLevels
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
+
+ end if
+
+
wdtz(nVertLevels+1) = 0.
+
do k=1,nVertLevels
tend_theta(k,iCell) = tend_theta(k,iCell)/areaCell(iCell) -rdzw(k)*(wdtz(k+1)-wdtz(k))
tend_theta(k,iCell) = tend_theta(k,iCell) + rho(k,iCell)*rt_diabatic(k,iCell)
@@ -2577,7 +2851,7 @@
!
if ( v_theta_eddy_visc2 > 0.0 ) then
- if (mix_full) then
+ if (config_mix_full) then
do iCell = 1, grid % nCellsSolve
do k=2,nVertLevels-1
@@ -3081,6 +3355,6 @@
end do
end do
- end subroutine kessler
+ end subroutine kessler
end module time_integration
Modified: branches/atmos_physics/src/core_nhyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/mpas_interface.F 2010-08-03 18:36:36 UTC (rev 457)
+++ branches/atmos_physics/src/core_nhyd_atmos/mpas_interface.F 2010-08-03 18:39:56 UTC (rev 458)
@@ -34,6 +34,7 @@
call init_coupled_diagnostics( block % time_levs(1) % state, mesh)
call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh) ! ok for nonhydrostatic model
call initialize_advection_rk(mesh)
+ call initialize_deformation_weights(mesh)
#ifdef DO_PHYSICS
!check that all the physics options are correctly defined and that at least one physics
@@ -84,9 +85,7 @@
#ifdef DO_PHYSICS
call physics_timetracker(itimestep)
if(l_physics) call physics_driver(domain,itimestep)
- write(6,*) '--- back to mpas_interface:'
call timestep(domain, dt, itimestep)
- write(6,*) '--- end timestep:'
write(6,*)
#else
call timestep(domain, dt)
</font>
</pre>