<p><b>laura@ucar.edu</b> 2010-05-20 14:20:35 -0600 (Thu, 20 May 2010)</p><p>Do not remove the f90 files when make clean.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_hyd_atmos/Makefile
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/Makefile 2010-05-20 20:04:14 UTC (rev 292)
+++ branches/atmos_physics/src/core_hyd_atmos/Makefile 2010-05-20 20:20:35 UTC (rev 293)
@@ -24,5 +24,5 @@
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
- $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
- $(RM) $*.f90
+ $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators -I../core_hyd_phys
+# $(RM) $*.f90
Modified: branches/atmos_physics/src/core_hyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/Registry 2010-05-20 20:04:14 UTC (rev 292)
+++ branches/atmos_physics/src/core_hyd_atmos/Registry 2010-05-20 20:20:35 UTC (rev 293)
@@ -24,7 +24,8 @@
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
-namelist real physics config_dt_physics 1800.
+namelist integer physics config_n_physics 01
+namelist integer physics config_n_microp 01
#
# dim type name_in_file name_in_code
@@ -108,11 +109,17 @@
var real qv ( nVertLevels nCells Time ) iro qv scalars moist
var real qc ( nVertLevels nCells Time ) iro qc scalars moist
var real qr ( nVertLevels nCells Time ) iro qr scalars moist
+var real qi ( nVertLevels nCells Time ) iro qi scalars moist
+var real qs ( nVertLevels nCells Time ) iro qs scalars moist
+var real qg ( nVertLevels nCells Time ) iro qg scalars moist
+var real qnr ( nVertLevels nCells Time ) iro qnr scalars number
+var real qni ( nVertLevels nCells Time ) iro qni scalars number
#var real tracers ( nTracers nVertLevels nCells Time ) iro tracers - -
# state variables diagnosed from prognostic state
var real h ( nVertLevels nCells Time ) ro h - -
var real ww ( nVertLevelsP1 nCells Time ) ro ww - -
+var real w ( nVertLevelsP1 nCells Time ) ro w - -
var real pressure ( nVertLevelsP1 nCells Time ) ro pressure - -
var real geopotential ( nVertLevelsP1 nCells Time ) ro geopotential - -
var real alpha ( nVertLevels nCells Time ) iro alpha - -
@@ -152,6 +159,11 @@
var real qv_old ( nVertLevels nCells ) - qv_old scalars_old moist_old
var real qc_old ( nVertLevels nCells ) - qc_old scalars_old moist_old
var real qr_old ( nVertLevels nCells ) - qr_old scalars_old moist_old
+var real qi_old ( nVertLevels nCells ) - qi_old scalars_old moist_old
+var real qs_old ( nVertLevels nCells ) - qs_old scalars_old moist_old
+var real qg_old ( nVertLevels nCells ) - qg_old scalars_old moist_old
+var real qnr_old ( nVertLevels nCells Time ) - qnr_old scalars_old number_old
+var real qni_old ( nVertLevels nCells Time ) - qni_old scalars_old number_old
#var real tracers_old ( nTracers nVertLevels nCells ) - tracers_old - -
# Space needed for advection
@@ -161,3 +173,35 @@
# Arrays required for reconstruction of velocity field
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
+#==================================================================================================
+#... PARAMETERIZATION OF CONVECTION:
+
+# cubot :
+# cutop :
+# pratec : convective precipitation rate (mm s-1)
+# raincv : time-step convective precipitation (mm)
+# rthcuten: coupled-tendency of potential temperature due to cumulus convection (Pa K s-1)
+# rqvcuten: coupled-tendency of water vapor mixing ratio due to cumulus convection (Pa kg/kg s-1)
+# rqccuten: coupled-tendency of cloud water mixing ratio due to cumulus convection (Pa kg/kg s-1)
+# rqrcuten: coupled-tendency of rain mixing ratio due to cumulus convection (Pa kg/kg s-1)
+# rqicuten: coupled-tendency of cloud ice mixing ratio due to cumulus convection (Pa kg/kg s-1)
+# rqscuten: coupled-tendency of snow mixing ratio due to cumulus convection (Pa kg/kg s-1)
+# wavg : average vertical velocity (KF scheme only) (m s-1)
+
+var real cubot ( nCells Time ) r cubot - -
+var real cutop ( nCells Time ) r cutop - -
+
+# PRECIPITATION:
+var real pratec ( nCells Time ) r pratec - -
+var real raincv ( nCells Time ) r raincv - -
+
+# TENDENCIES:
+var real rthcuten ( nVertLevels nCells Time ) iro rthcuten - -
+var real rqvcuten ( nVertLevels nCells Time ) iro rqvcuten - -
+var real rqccuten ( nVertLevels nCells Time ) iro rqccuten - -
+var real rqrcuten ( nVertLevels nCells Time ) iro rqrcuten - -
+var real rqicuten ( nVertLevels nCells Time ) iro rqicuten - -
+var real rqscuten ( nVertLevels nCells Time ) iro rqscuten - -
+var real wavg ( nVertLevels nCells Time ) iro wavg - -
+
+#==================================================================================================
Modified: branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F 2010-05-20 20:04:14 UTC (rev 292)
+++ branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F 2010-05-20 20:20:35 UTC (rev 293)
@@ -6,11 +6,19 @@
use dmpar
use vector_reconstruction
+#ifdef DO_PHYSICS
+ use module_microphysics_driver
+#endif
contains
+#ifdef DO_PHYSICS
+ subroutine timestep(domain, dt, itimestep)
+#else
subroutine timestep(domain, dt)
+#endif
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step
!
@@ -27,8 +35,16 @@
type (block_type), pointer :: block
+#ifdef DO_PHYSICS
+ integer, intent(in):: itimestep
+#endif
+
if (trim(config_time_integration) == 'SRK3') then
+#ifdef DO_PHYSICS
+ call srk3(domain, dt, itimestep)
+#else
call srk3(domain, dt)
+#endif
else
write(0,*) 'Unknown time integration option '//trim(config_time_integration)
write(0,*) 'Currently, only ''SRK3'' is supported.'
@@ -40,11 +56,16 @@
block % time_levs(2) % state % xtime % scalar = block % time_levs(1) % state % xtime % scalar + dt
block => block % next
end do
+ write(0,*) '--- end subroutine timestep:'
end subroutine timestep
+#ifdef DO_PHYSICS
+ subroutine srk3(domain, dt, itimestep)
+#else
subroutine srk3(domain, dt)
+#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step using
! time-split RK3 scheme
@@ -77,6 +98,10 @@
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
+#ifdef DO_PHYSICS
+ integer, intent(in):: itimestep
+#endif
+
!
! Initialize time_levs(2) with state at current time
! Initialize RK weights
@@ -315,6 +340,14 @@
! END RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Compute vertical velocity diagnostic
+ !
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_w( block % time_levs(2) % state, block % time_levs(1) % state, block % mesh, dt )
+ block => block % next
+ end do
!
! Compute full velocity vectors at cell centers
@@ -361,6 +394,32 @@
end if
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! LDF begin (04-26-2010):
+ ! CALL TO CLOUD MICROPHYISCS SCHEMES STARTS HERE:
+ ! LDF end.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!#ifdef DO_PHYSICS
+! if(config_mp_physics /=0) then
+
+! block => domain % blocklist
+! do while(associated(block))
+! call microphysics_driver ( block % mesh, block % time_levs(2) % state, itimestep)
+! block => block % next
+! end do
+
+! endif
+!#endif
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! LDF begin (04-26-2010):
+ ! CALL TO CLOUD MICROPHYISCS SCHEMES ENDS HERE:
+ ! LDF end.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ write(0,*) '--- end subroutine srk3:'
+
end subroutine srk3
!------------------------------------------------------------------------------------------------------------------
@@ -683,10 +742,10 @@
if ( h_mom_eddy_visc4 > 0.0 ) then
- allocate(delsq_divergence(nVertLevels, nCells))
- allocate(delsq_u(nVertLevels, nEdges))
- allocate(delsq_circulation(nVertLevels, nVertices))
- allocate(delsq_vorticity(nVertLevels, nVertices))
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
delsq_u(:,:) = 0.0
@@ -861,7 +920,7 @@
if ( h_theta_eddy_visc4 > 0.0 ) then
- allocate(delsq_theta(nVertLevels, nCells))
+ allocate(delsq_theta(nVertLevels, nCells+1))
delsq_theta(:,:) = 0.
@@ -1143,7 +1202,7 @@
p0 = 1.e+05
ptop = pressure(nVertLevels+1,1)
smext = 0.1
- smdiv = 0.1
+ smdiv = 0.0
! write(0,*) ' ptop in advance_dynamics ',ptop
@@ -1531,9 +1590,9 @@
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension( num_scalars, grid % nEdges) :: h_flux
- real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: v_flux, v_flux_upwind, s_update
- real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: scale_out, scale_in
+ real (kind=RKIND), dimension( num_scalars, grid % nEdges+1) :: h_flux
+ real (kind=RKIND), dimension( num_scalars, grid % nCells+1, 2 ) :: v_flux, v_flux_upwind, s_update
+ real (kind=RKIND), dimension( num_scalars, grid % nCells+1, 2 ) :: scale_out, scale_in
real (kind=RKIND), dimension( num_scalars ) :: s_max, s_min, s_max_update, s_min_update
integer :: nVertLevels, km0, km1, ktmp, kcp1, kcm1
@@ -2114,4 +2173,84 @@
end subroutine compute_solve_diagnostics
+
+ subroutine compute_w (s_new, s_old, grid, dt )
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute (diagnose) vertical velocity (used by physics)
+ !
+ ! Input: s_new - current model state
+ ! s_old - previous model state
+ ! grid - grid metadata
+ ! dt - timestep between new and old
+ !
+ ! Output: w (m/s)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_state), intent(inout) :: s_new
+ type (grid_state), intent(in) :: s_old
+ type (grid_meta), intent(inout) :: grid
+
+ real (kind=RKIND), intent(in) :: dt
+
+ real (kind=RKIND), dimension(:,:), pointer :: geo_new, geo_old, u, ww, h, w
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, rdnw, fnm, fnp
+
+ integer :: iEdge, iCell, k, cell1, cell2
+ real (kind=RKIND), dimension( grid % nVertlevels + 1 ) :: wdwn
+ real (kind=RKIND) :: flux
+
+ geo_new => s_new % geopotential % array
+ geo_old => s_old % geopotential % array
+ u => s_new % u % array
+ ww => s_new % ww % array
+ h => s_new % h % array
+ w => s_new % w % array
+ dvEdge => grid % dvEdge % array
+ rdnw => grid % rdnw % array
+ fnm => grid % fnm % array
+ fnp => grid % fnp % array
+
+ w = 0.
+
+ do iCell=1, grid % nCellsSolve
+ wdwn(1) = 0.
+ do k=2,grid % nVertlevels + 1
+ wdwn(k) = (0.5*(ww(k,iCell)+ww(k-1,iCell))/h(k-1,iCell)) &
+ *rdnw(k-1)*(geo_new(k,iCell)-geo_new(k-1,iCell))
+ enddo
+ w(1,iCell) = 0.
+ do k=2, grid % nVertLevels
+ w(k,iCell) = (((geo_new(k,iCell)-geo_old(k,iCell))/dt)+ &
+ fnm(k)*wdwn(k+1)+fnp(k)*wdwn(k))/gravity
+ enddo
+ k = grid % nVertLevels + 1
+ w(k,iCell) = ((geo_new(k,iCell)-geo_old(k,iCell))/dt)/gravity
+ enddo
+
+ do iEdge=1, grid % nEdges
+ cell1 = grid % cellsOnEdge % array (1,iEdge)
+ cell2 = grid % cellsOnEdge % array (2,iEdge)
+ if (cell1 <= grid % nCellsSolve .or. cell2 <= grid % nCellsSolve ) then
+ do k=2, grid % nVertLevels
+ flux = 0.25*(u(k,iEdge)+u(k-1,iEdge)*(geo_new(k,cell2)-geo_new(k,cell1))) * dvEdge(iEdge) / gravity
+ w(k,cell1) = w(k,cell1) + flux/ grid % AreaCell % array (cell1)
+ w(k,cell2) = w(k,cell2) + flux/ grid % AreaCell % array (cell2)
+ enddo
+ k = 1
+ flux = 0.5*(1.5*u(1,iEdge)-0.5*u(2,iEdge)*(geo_new(k,cell2)-geo_new(k,cell1))) * dvEdge(iEdge) / gravity
+ w(k,cell1) = w(k,cell1) + flux/ grid % AreaCell % array (cell1)
+ w(k,cell2) = w(k,cell2) + flux/ grid % AreaCell % array (cell2)
+
+ k = grid % nVertLevels + 1
+ flux = 0.5*(1.5*u(k-1,iEdge)-0.5*u(k-2,iEdge)*(geo_new(k,cell2)-geo_new(k,cell1))) * dvEdge(iEdge) / gravity
+ w(k,cell1) = w(k,cell1) + flux/ grid % AreaCell % array (cell1)
+ w(k,cell2) = w(k,cell2) + flux/ grid % AreaCell % array (cell2)
+
+ end if
+ end do
+
+ end subroutine compute_w
+
end module time_integration
Modified: branches/atmos_physics/src/core_hyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/mpas_interface.F 2010-05-20 20:04:14 UTC (rev 292)
+++ branches/atmos_physics/src/core_hyd_atmos/mpas_interface.F 2010-05-20 20:20:35 UTC (rev 293)
@@ -18,6 +18,10 @@
use grid_types
use advection
use time_integration
+#ifdef DO_PHYSICS
+ use module_physics_init
+ use module_physics_manager
+#endif
implicit none
@@ -32,6 +36,11 @@
call init_reconstruct(mesh)
call reconstruct(block % time_levs(1) % state, mesh)
+#ifdef DO_PHYSICS
+ call physics_wrf_allocate(mesh)
+! call physics_init(mesh, block % time_levs(1) % state)
+#endif
+
end subroutine mpas_init
@@ -49,24 +58,49 @@
end subroutine mpas_query
-subroutine mpas_timestep(domain, itimestep, dt)
+subroutine mpas_timestep(domain, itimestep, dt &
+#ifdef DO_PHYSICS
+ ,n_physics &
+#endif
+ )
use grid_types
use time_integration
+#ifdef DO_PHYSICS
+ use module_physics_driver
+ use module_physics_manager
+#endif
implicit none
type (domain_type), intent(inout) :: domain
integer, intent(in) :: itimestep
real (kind=RKIND), intent(in) :: dt
+#ifdef DO_PHYSICS
+ integer, intent(in) :: n_physics
+ logical:: l_physics
+#endif
+#ifdef DO_PHYSICS
+! l_physics = .false.
+! call physics_timetracker(itimestep,l_physics)
+! if(l_physics) call physics_driver(domain,itimestep)
+ call timestep(domain, dt, itimestep)
+#else
call timestep(domain, dt)
+#endif
end subroutine mpas_timestep
subroutine mpas_finalize()
+#ifdef DO_PHYSICS
+ use module_physics_manager
+ implicit none
+ call physics_wrf_deallocate
+#else
implicit none
+#endif
end subroutine mpas_finalize
</font>
</pre>