<p><b>duda</b> 2012-02-27 15:31:30 -0700 (Mon, 27 Feb 2012)</p><p>Merge changes from atmos_physics branch; only affects atmosphere cores.<br>
<br>
<br>
M namelist.input.nhyd_atmos_squall<br>
M src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F<br>
M src/core_init_nhyd_atmos/Registry<br>
M src/core_atmos_physics/mpas_atmphys_manager.F<br>
M src/core_atmos_physics/mpas_atmphys_driver_convection_deep.F<br>
M src/core_nhyd_atmos/mpas_atm_test_cases.F<br>
M src/core_nhyd_atmos/mpas_atm_time_integration.F<br>
M src/core_nhyd_atmos/Registry<br>
M namelist.input.nhyd_atmos<br>
M namelist.input.nhyd_atmos_mtn_wave<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/namelist.input.nhyd_atmos
===================================================================
--- trunk/mpas/namelist.input.nhyd_atmos 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/namelist.input.nhyd_atmos 2012-02-27 22:31:30 UTC (rev 1537)
@@ -1,5 +1,4 @@
&nhyd_model
- config_test_case = 0
config_time_integration = 'SRK3'
config_dt = 450.0
config_start_time = '2010-10-23_00:00:00'
Modified: trunk/mpas/namelist.input.nhyd_atmos_mtn_wave
===================================================================
--- trunk/mpas/namelist.input.nhyd_atmos_mtn_wave 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/namelist.input.nhyd_atmos_mtn_wave 2012-02-27 22:31:30 UTC (rev 1537)
@@ -1,49 +1,72 @@
&nhyd_model
- config_test_case = 6
config_time_integration = 'SRK3'
- config_dt = 6.
- config_ntimesteps = 3000
- config_output_interval = 100
+ config_dt = 6.0
+ config_start_time = '0000-01-01_00:00:00'
+ config_run_duration = '05:00:00'
config_number_of_sub_steps = 6
config_h_mom_eddy_visc2 = 10.
config_h_mom_eddy_visc4 = 0.
- config_v_mom_eddy_visc2 = 10.
+ config_v_mom_eddy_visc2 = 10.0
config_h_theta_eddy_visc2 = 10.
config_h_theta_eddy_visc4 = 0.
config_v_theta_eddy_visc2 = 10.
+ config_horiz_mixing = '2d_fixed'
config_theta_adv_order = 3
config_w_adv_order = 3
+ config_u_vadv_order = 3
+ config_w_vadv_order = 3
+ config_theta_vadv_order = 3
+ config_coef_3rd_order = 0.25
config_scalar_advection = .false.
- config_positive_definite = .false.
- config_monotonic = .false.
- config_mix_full = .false.
- config_horiz_mixing = '2d_fixed'
- config_coef_3rd_order = 0.25
- config_epssm = 0.2
+ config_epssm = 0.1
+ config_smdiv = 0.1
config_newpx = .false.
/
-
&damping
- config_zd = 22000.0
- config_xnutr = 0.0
+ config_zd = 10500.0
+ config_xnutr = 0.1
/
&dimensions
- config_nvertlevels = 26
+ config_nvertlevels = 70
/
&io
- config_input_name = 'grid.nc'
+ config_input_name = 'mtn_wave_init.nc'
config_output_name = 'output.nc'
config_restart_name = 'restart.nc'
+ config_output_interval = '00:30:00'
+ config_frames_per_outfile = 0
/
&restart
- config_restart_interval = 3000
+ config_restart_interval = '1_00:00:00'
config_do_restart = .false.
- config_restart_time = 1036800.0
/
&physics
+ config_frac_seaice = .false.
+ config_sfc_albedo = .false.
+ config_sst_update = .false.
+ config_sstdiurn_update = .false.
+ config_deepsoiltemp_update = .false.
+
+ config_n_microp = 5
+
+ config_radtlw_interval = '00:00:00'
+ config_radtsw_interval = '00:00:00'
+ config_conv_interval = 'none'
+ config_pbl_interval = 'none'
+
+ config_microp_scheme = 'off'
+ config_conv_shallow_scheme = 'off'
+ config_conv_deep_scheme = 'off'
+ config_eddy_scheme = 'off'
+ config_lsm_scheme = 'off'
+ config_pbl_scheme = 'off'
+ config_radt_cld_scheme = 'off'
+ config_radt_lw_scheme = 'off'
+ config_radt_sw_scheme = 'off'
+ config_sfclayer_scheme = 'off'
/
Modified: trunk/mpas/namelist.input.nhyd_atmos_squall
===================================================================
--- trunk/mpas/namelist.input.nhyd_atmos_squall 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/namelist.input.nhyd_atmos_squall 2012-02-27 22:31:30 UTC (rev 1537)
@@ -1,43 +1,74 @@
&nhyd_model
- config_test_case = 4
config_time_integration = 'SRK3'
- config_dt = 6.
- config_ntimesteps = 600
- config_output_interval = 100
+ config_dt = 3.0
+ config_start_time = '0000-01-01_00:00:00'
+ config_run_duration = '02:00:00'
config_number_of_sub_steps = 6
config_h_mom_eddy_visc2 = 500.
config_h_mom_eddy_visc4 = 0.
config_v_mom_eddy_visc2 = 500.0
config_h_theta_eddy_visc2 = 500.
- config_h_theta_eddy_visc4 = 00.
- config_v_theta_eddy_visc2 = 500.0
- config_theta_adv_order = 2
- config_scalar_adv_order = 2
- config_positive_definite = .false.
- config_monotonic = .false.
- config_newpx = .false.
+ config_h_theta_eddy_visc4 = 0.
+ config_v_theta_eddy_visc2 = 500.
+ config_horiz_mixing = '2d_fixed'
+ config_theta_adv_order = 3
+ config_w_adv_order = 3
+ config_u_vadv_order = 3
+ config_w_vadv_order = 3
+ config_theta_vadv_order = 3
+ config_coef_3rd_order = 0.25
+ config_epssm = 0.1
+ config_smdiv = 0.1
+ config_mix_full = .false.
+ config_monotonic = .true.
/
&damping
- config_zd = 22000.0
+ config_zd = 20000.0
config_xnutr = 0.0
/
&dimensions
- config_nvertlevels = 26
+ config_nvertlevels = 40
/
&io
- config_input_name = 'grid.nc'
+ config_input_name = 'supercell_init.nc'
config_output_name = 'output.nc'
config_restart_name = 'restart.nc'
+ config_output_interval = '00:30:00'
+ config_frames_per_outfile = 0
/
&restart
- config_restart_interval = 3000
+ config_restart_interval = '1_00:00:00'
config_do_restart = .false.
- config_restart_time = 1036800.0
/
&physics
+ config_frac_seaice = .false.
+ config_sfc_albedo = .false.
+ config_sst_update = .false.
+ config_sstdiurn_update = .false.
+ config_deepsoiltemp_update = .false.
+
+ config_n_microp = 1
+
+ config_radtlw_interval = '00:30:00'
+ config_radtsw_interval = '00:30:00'
+ config_conv_interval = 'none'
+ config_pbl_interval = 'none'
+
+ config_microp_scheme = 'kessler'
+ config_conv_shallow_scheme = 'off'
+ config_conv_deep_scheme = 'off'
+ config_eddy_scheme = 'off'
+ config_lsm_scheme = 'off'
+ config_pbl_scheme = 'off'
+ config_radt_cld_scheme = 'off'
+ config_radt_lw_scheme = 'off'
+ config_radt_sw_scheme = 'off'
+ config_sfclayer_scheme = 'off'
/
+
+
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_convection_deep.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_convection_deep.F 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_convection_deep.F 2012-02-27 22:31:30 UTC (rev 1537)
@@ -196,7 +196,7 @@
write(0,*) '--- enter subroutine kf_eta_cps:'
call kf_eta_cps ( &
dt = dt_dyn , ktau = itimestep , &
- areaCell = area_p , cudt = dt_cu , &
+ areaCell = area_p , cudt = cudt , &
curr_secs = curr_secs , adapt_step_flag = adapt_step_flag , &
rho = rho_p , raincv = raincv_p , &
pratec = pratec_p , nca = nca_p , &
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_manager.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_manager.F 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_manager.F 2012-02-27 22:31:30 UTC (rev 1537)
@@ -290,7 +290,7 @@
if(ierr /= 0) &
call physics_error_fatal('subroutine physics_run_init: error defining dt_pbl')
- elseif(trim(config_conv_interval) == "none") then
+ elseif(trim(config_pbl_interval) == "none") then
dt_pbl = config_dt
else
Modified: trunk/mpas/src/core_init_nhyd_atmos/Registry
===================================================================
--- trunk/mpas/src/core_init_nhyd_atmos/Registry 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/src/core_init_nhyd_atmos/Registry 2012-02-27 22:31:30 UTC (rev 1537)
@@ -1,24 +1,24 @@
%
% namelist type namelist_record name default_value
%
-namelist integer nhyd_model config_test_case 5
+namelist integer nhyd_model config_test_case 7
namelist integer nhyd_model config_calendar_type MPAS_GREGORIAN
namelist character nhyd_model config_start_time none
namelist character nhyd_model config_stop_time none
-namelist integer nhyd_model config_theta_adv_order 2
-namelist real nhyd_model config_coef_3rd_order 1.0
+namelist integer nhyd_model config_theta_adv_order 3
+namelist real nhyd_model config_coef_3rd_order 0.25
namelist integer dimensions config_nvertlevels 26
namelist integer dimensions config_nsoillevels 4
namelist integer dimensions config_nfglevels 27
namelist integer dimensions config_nfgsoillevels 4
namelist integer dimensions config_months 12
-namelist character data_sources config_geog_data_path /data3/mp/wrfhelp/WPS_GEOG/
+namelist character data_sources config_geog_data_path /mmm/users/wrfhelp/WPS_GEOG/
namelist character data_sources config_met_prefix FILE
namelist character data_sources config_sfc_prefix FILE
namelist integer data_sources config_fg_interval 21600
-namelist real vertical_grid config_ztop 24000.0
+namelist real vertical_grid config_ztop 28000.0
namelist integer vertical_grid config_nsmterrain 2
-namelist logical vertical_grid config_smooth_surfaces true
+namelist logical vertical_grid config_smooth_surfaces false
namelist logical preproc_stages config_static_interp true
namelist logical preproc_stages config_vertical_grid true
namelist logical preproc_stages config_met_interp true
@@ -26,7 +26,7 @@
namelist logical preproc_stages config_frac_seaice false
namelist character io config_input_name grid.nc
namelist character io config_sfc_update_name sfc_update.nc
-namelist character io config_output_name output.nc
+namelist character io config_output_name init.nc
namelist character io config_restart_name restart.nc
namelist character io config_decomp_file_prefix graph.info.part.
namelist integer io config_frames_per_outfile 0
Modified: trunk/mpas/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- trunk/mpas/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-02-27 22:31:30 UTC (rev 1537)
@@ -35,6 +35,28 @@
type (block_type), pointer :: block_ptr
+
+ !
+ ! Do some quick checks to make sure compile options are compatible with the chosen test case
+ !
+ if (config_test_case == 6) then
+#ifndef ROTATED_GRID
+ write(0,*) '*** ERROR ***'
+ write(0,*) 'To initialize and run the mountain wave test case (case 6),'
+ write(0,*) ' you must compile with -DROTATE_GRID added to the specification'
+ write(0,*) ' of MODEL_FORMULATION at the top of the Makefile.'
+ call mpas_dmpar_abort(domain % dminfo)
+#endif
+ else
+#ifdef ROTATED_GRID
+ write(0,*) '*** ERROR ***'
+ write(0,*) 'Only test case 6 should use code compiled with -DROTATE_GRID specified in the Makefile.'
+ call mpas_dmpar_abort(domain % dminfo)
+#endif
+ end if
+
+
+
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 '
@@ -83,7 +105,7 @@
do while (associated(block_ptr))
call init_atm_test_case_gfs(domain % dminfo, block_ptr % mesh, block_ptr % fg, block_ptr % state % time_levs(1) % state, &
block_ptr % diag, config_test_case, block_ptr % parinfo)
- call physics_initialize_real(block_ptr % mesh, block_ptr % fg)
+ if (config_met_interp) call physics_initialize_real(block_ptr % mesh, block_ptr % fg)
block_ptr => block_ptr % next
end do
@@ -1877,7 +1899,11 @@
+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)
+#ifdef ROTATED_GRID
+ u(k,i) = sin(grid % angleEdge % array(i)) * (u(k,i) - us)
+#else
+ u(k,i) = cos(grid % angleEdge % array(i)) * (u(k,i) - us)
+#endif
end do
end if
end do
Modified: trunk/mpas/src/core_nhyd_atmos/Registry
===================================================================
--- trunk/mpas/src/core_nhyd_atmos/Registry 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/src/core_nhyd_atmos/Registry 2012-02-27 22:31:30 UTC (rev 1537)
@@ -1,9 +1,9 @@
%
% namelist type namelist_record name default_value
%
-namelist integer nhyd_model config_test_case 5
+namelist integer nhyd_model config_test_case 0
namelist character nhyd_model config_time_integration SRK3
-namelist real nhyd_model config_dt 172.8
+namelist real nhyd_model config_dt 600.0
namelist integer nhyd_model config_calendar_type MPAS_GREGORIAN
namelist character nhyd_model config_start_time 0000-01-01_00:00:00
namelist character nhyd_model config_stop_time none
@@ -17,19 +17,19 @@
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 integer nhyd_model config_w_adv_order 3
+namelist integer nhyd_model config_theta_adv_order 3
+namelist integer nhyd_model config_scalar_adv_order 3
+namelist integer nhyd_model config_u_vadv_order 3
+namelist integer nhyd_model config_w_vadv_order 3
+namelist integer nhyd_model config_theta_vadv_order 3
+namelist integer nhyd_model config_scalar_vadv_order 3
+namelist real nhyd_model config_coef_3rd_order 0.25
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 real nhyd_model config_len_disp 120000.0
namelist real nhyd_model config_epssm 0.1
namelist real nhyd_model config_smdiv 0.1
namelist logical nhyd_model config_newpx false
@@ -38,7 +38,7 @@
namelist real damping config_zd 22000.0
namelist real damping config_xnutr 0.0
namelist integer dimensions config_nvertlevels 26
-namelist character io config_input_name grid.nc
+namelist character io config_input_name init.nc
namelist character io config_sfc_update_name sfc_update.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -333,14 +333,14 @@
namelist logical physics config_sstdiurn_update false
namelist logical physics config_deepsoiltemp_update false
-namelist integer physics config_n_physics 01
-namelist integer physics config_n_microp 01
-namelist integer physics config_n_conv 01
-namelist integer physics config_n_pbl 01
-namelist integer physics config_n_lsm 01
-namelist integer physics config_n_eddy 01
-namelist integer physics config_n_radt_lw 01
-namelist integer physics config_n_radt_sw 01
+namelist integer physics config_n_physics 1
+namelist integer physics config_n_microp 1
+namelist integer physics config_n_conv 1
+namelist integer physics config_n_pbl 1
+namelist integer physics config_n_lsm 1
+namelist integer physics config_n_eddy 1
+namelist integer physics config_n_radt_lw 1
+namelist integer physics config_n_radt_sw 1
namelist character physics config_radtlw_interval none
namelist character physics config_radtsw_interval none
Modified: trunk/mpas/src/core_nhyd_atmos/mpas_atm_test_cases.F
===================================================================
--- trunk/mpas/src/core_nhyd_atmos/mpas_atm_test_cases.F 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/src/core_nhyd_atmos/mpas_atm_test_cases.F 2012-02-27 22:31:30 UTC (rev 1537)
@@ -1865,7 +1865,11 @@
+zgrid(k,cell2)+zgrid(k+1,cell2))
u(k,i) = um
if(i == 1 ) grid % u_init % array(k) = u(k,i) - us
+#ifdef ROTATED_GRID
+ u(k,i) = sin(grid % angleEdge % array(i)) * (u(k,i) - us)
+#else
u(k,i) = cos(grid % angleEdge % array(i)) * (u(k,i) - us)
+#endif
end do
end if
end do
Modified: trunk/mpas/src/core_nhyd_atmos/mpas_atm_time_integration.F
===================================================================
--- trunk/mpas/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-02-27 21:54:55 UTC (rev 1536)
+++ trunk/mpas/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-02-27 22:31:30 UTC (rev 1537)
@@ -93,8 +93,6 @@
logical, parameter :: debug = .false.
! logical, parameter :: debug = .true.
logical, parameter :: debug_mass_conservation = .true.
-! logical, parameter :: do_microphysics = .true.
- logical, parameter :: do_microphysics = .false.
integer :: index_qc
real (kind=RKIND) :: domain_mass, scalar_mass, scalar_min, scalar_max
@@ -415,13 +413,6 @@
end do
#endif
-! if(do_microphysics) then
-! block => domain % blocklist
-! do while (associated(block))
-! call atm_qd_kessler( block % state % time_levs(1) % state, block % state % time_levs(2) % state, block % mesh, dt )
-! block => block % next
-! end do
-! end if
! if(debug) then
101 format(' local min, max scalar',i4,2(1x,e17.10))
@@ -2370,7 +2361,11 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1,nVertLevels
+#ifdef ROTATED_GRID
+ u_mix(k) = u(k,iEdge) - grid % u_init % array(k) * sin( grid % angleEdge % array(iEdge) )
+#else
u_mix(k) = u(k,iEdge) - grid % u_init % array(k) * cos( grid % angleEdge % array(iEdge) )
+#endif
end do
do k=2,nVertLevels-1
@@ -3449,161 +3444,5 @@
end subroutine atm_init_coupled_diagnostics
-! ------------------------
- subroutine atm_qd_kessler( state_old, state_new, diag, tend, grid, dt )
-
- implicit none
-
- type (state_type), intent(inout) :: state_old, state_new
- type (diag_type), intent(inout) :: diag
- type (tend_type), intent(inout) :: tend
- type (mesh_type), intent(inout) :: grid
- real (kind=RKIND), intent(in) :: dt
-
- real (kind=RKIND), dimension( grid % nVertLevels ) :: t, rho, p, dzu, qv, qc, qr, qc1, qr1
-
- integer :: k,iEdge,i,iCell,nz1
- real (kind=RKIND) :: p0,rcv
-
-
- write(0,*) ' in qd_kessler '
-
- p0 = 1.e+05
- rcv = rgas/(cp-rgas)
- nz1 = grid % nVertLevels
-
- do iCell = 1, grid % nCellsSolve
-
- do k = 1, grid % nVertLevels
-
- tend % rt_diabatic_tend % array(k,iCell) = state_new % theta_m % array(k,iCell)
-
- t(k) = state_new % theta_m % array(k,iCell)/(1. + 1.61*state_new % scalars % array(state_new % index_qv,k,iCell))
- rho(k) = grid % zz % array(k,iCell)*state_new % rho_zz % array(k,iCell)
- p(k) = diag % exner % array(k,iCell)
- qv(k) = max(0.0_RKIND,state_new % scalars % array(state_new % index_qv,k,iCell))
- qc(k) = max(0.0_RKIND,state_new % scalars % array(state_new % index_qc,k,iCell))
- qr(k) = max(0.0_RKIND,state_new % scalars % array(state_new % index_qr,k,iCell))
- qc1(k) = max(0.0_RKIND,state_old % scalars % array(state_old % index_qc,k,iCell))
- qr1(k) = max(0.0_RKIND,state_old % scalars % array(state_old % index_qr,k,iCell))
- dzu(k) = grid % dzu % array(k)
-
- end do
-
- call atm_kessler( t,qv,qc,qc1,qr,qr1,rho,p,dt,dzu,nz1, 1)
-
- do k = 1, grid % nVertLevels
-
- state_new % theta_m % array(k,iCell) = t(k)*(1.+1.61*qv(k))
- tend % rt_diabatic_tend % array(k,iCell) = state_new % rho_zz % array(k,iCell) * &
- (state_new % theta_m % array(k,iCell) - tend % rt_diabatic_tend % array(k,iCell))/dt
- diag % rtheta_p % array(k,iCell) = state_new % rho_zz % array(k,iCell) * state_new % theta_m % array(k,iCell) &
- - diag % rtheta_base % array(k,iCell)
- state_new % scalars % array(state_new % index_qv,k,iCell) = qv(k)
- state_new % scalars % array(state_new % index_qc,k,iCell) = qc(k)
- state_new % scalars % array(state_new % index_qr,k,iCell) = qr(k)
-
- diag % exner % array(k,iCell) = &
- ( grid % zz % array(k,iCell)*(rgas/p0) * ( &
- diag % rtheta_p % array(k,iCell) &
- + diag % rtheta_base % array(k,iCell) ) )**rcv
-
- diag % pressure_p % array(k,iCell) = &
- grid % zz % array(k,iCell) * rgas * ( &
- diag % exner % array(k,iCell)*diag % rtheta_p % array(k,iCell) &
- +diag % rtheta_base % array(k,iCell) * &
- (diag % exner % array(k,iCell) - diag % exner_base % array(k,iCell)) )
- end do
-
- end do
-
- write(0,*) ' exiting qd_kessler '
-
- end subroutine atm_qd_kessler
-
-!-----------------------------------------------------------------------
- subroutine atm_kessler( t1t, qv1t, qc1t, qc1, qr1t, qr1, &
- rho, pii, dt, dzu, nz1, nx )
-!-----------------------------------------------------------------------
-!
- implicit none
- integer :: nx, nz1
- real (kind=RKIND) :: t1t (nz1,nx), qv1t(nz1,nx), qc1t(nz1,nx), &
- qr1t(nz1,nx), qc1 (nz1,nx), qr1 (nz1,nx), &
- rho (nz1,nx), pii (nz1,nx), dzu(nz1)
- integer, parameter :: mz=200
- real (kind=RKIND) :: qrprod(mz), prod (mz), rcgs( mz), rcgsi (mz) &
- ,ern (mz), vt (mz), vtden(mz), gam (mz) &
- ,r (mz), rhalf(mz), velqr(mz), buoycy(mz) &
- ,pk (mz), pc (mz), f0 (mz), qvs (mz)
-
- real (kind=RKIND) :: c1, c2, c3, c4, f5, mxfall, dtfall, fudge, dt, velu, veld, artemp, artot
- real (kind=RKIND) :: cp, product, ackess, ckess, fvel, f2x, xk, xki, psl
- integer :: nfall
- integer :: i,k,n
-
- ackess = 0.001
- ckess = 2.2
- fvel = 36.34
- f2x = 17.27
- f5 = 237.3*f2x*2.5e6/1003.
- xk = .2875
- xki = 1./xk
- psl = 1000.
-
- do k=1,nz1
- r(k) = 0.001*rho(k,1)
- rhalf(k) = sqrt(rho(1,1)/rho(k,1))
- pk(k) = pii(k,1)
- pc(k) = 3.8/(pk(k)**xki*psl)
- f0(k) = 2.5e6/(1003.*pk(k))
- end do
-!
- do i=1,nx
- do k=1,nz1
- qrprod(k) = qc1t(k,i) &
- -(qc1t(k,i)-dt*max(ackess*(qc1(k,i)-.001), &
- 0.0_RKIND))/(1.+dt*ckess*qr1(k,i)**.875)
- velqr(k) = (qr1(k,i)*r(k))**1.1364*rhalf(k)
- qvs(k) = pc(k)*exp(f2x*(pk(k)*t1t(k,i)-273.) &
- /(pk(k)*t1t(k,i)- 36.))
- end do
- velu = (qr1(2,i)*r(2))**1.1364*rhalf(2)
- veld = (qr1(1,i)*r(1))**1.1364*rhalf(1)
- qr1t(1,i) = qr1t(1,i)+dt*(velu-veld)*fvel/(r(1)*dzu(2))
- do k=2,nz1-1
- qr1t(k,i) = qr1t(k,i)+dt*fvel/r(k) &
- *.5*((velqr(k+1)-velqr(k ))/dzu(k+1) &
- +(velqr(k )-velqr(k-1))/dzu(k ))
- end do
- qr1t(nz1,i) = qr1t(nz1,i)-dt*fvel*velqr(nz1-1) &
- /(r(nz1)*dzu(nz1)*(1.+1.))
- artemp = 36340.*(.5*(velqr(2)+velqr(1))+veld-velu)
- artot = artot+dt*artemp
- do k=1,nz1
- qc1t(k,i) = max(qc1t(k,i)-qrprod(k),0.0_RKIND)
- qr1t(k,i) = max(qr1t(k,i)+qrprod(k),0.0_RKIND)
- prod(k) = (qv1t(k,i)-qvs(k))/(1.+qvs(k)*f5 &
- /(pk(k)*t1t(k,i)-36.)**2)
- end do
- do k=1,nz1
- ern(k) = min(dt*(((1.6+124.9*(r(k)*qr1t(k,i))**.2046) &
- *(r(k)*qr1t(k,i))**.525)/(2.55e6*pc(k) &
- /(3.8 *qvs(k))+5.4e5))*(dim(qvs(k),qv1t(k,i)) &
- /(r(k)*qvs(k))), &
- max(-prod(k)-qc1t(k,i),0.0_RKIND),qr1t(k,i))
- end do
- do k=1,nz1
- buoycy(k) = f0(k)*(max(prod(k),-qc1t(k,i))-ern(k))
- qv1t(k,i) = max(qv1t(k,i) &
- -max(prod(k),-qc1t(k,i))+ern(k),0.0_RKIND)
- qc1t(k,i) = qc1t(k,i)+max(prod(k),-qc1t(k,i))
- qr1t(k,i) = qr1t(k,i)-ern(k)
- t1t (k,i) = t1t (k,i)+buoycy(k)
- end do
- end do
-
- end subroutine atm_kessler
-
end module atm_time_integration
</font>
</pre>