Home

Context Navigation

source: palm/trunk/SOURCE/microphysics.f90 @ 1684

Last change on this file since 1684 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 75.4 KB

Rev	Line
[1682]	1	!> @file microphysics.f90
[1093]	2	!--------------------------------------------------------------------------------!
	3	! This file is part of PALM.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the terms
	6	! of the GNU General Public License as published by the Free Software Foundation,
	7	! either version 3 of the License, or (at your option) any later version.
	8	!
	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	12	!
	13	! You should have received a copy of the GNU General Public License along with
	14	! PALM. If not, see <http://www.gnu.org/licenses/>.
	15	!
[1310]	16	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1093]	17	!--------------------------------------------------------------------------------!
	18	!
[1000]	19	! Current revisions:
[1092]	20	! ------------------
[1683]	21	!
	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: microphysics.f90 1683 2015-10-07 23:57:51Z knoop $
	26	!
[1683]	27	! 1682 2015-10-07 23:56:08Z knoop
	28	! Code annotations made doxygen readable
	29	!
[1647]	30	! 1646 2015-09-02 16:00:10Z hoffmann
	31	! Bugfix: Wrong computation of d_mean.
	32	!
[1362]	33	! 1361 2014-04-16 15:17:48Z hoffmann
	34	! Bugfix in sedimentation_rain: Index corrected.
	35	! Vectorized version of adjust_cloud added.
	36	! Little reformatting of the code.
	37	!
[1354]	38	! 1353 2014-04-08 15:21:23Z heinze
	39	! REAL constants provided with KIND-attribute
	40	!
[1347]	41	! 1346 2014-03-27 13:18:20Z heinze
	42	! Bugfix: REAL constants provided with KIND-attribute especially in call of
	43	! intrinsic function like MAX, MIN, SIGN
	44	!
[1335]	45	! 1334 2014-03-25 12:21:40Z heinze
	46	! Bugfix: REAL constants provided with KIND-attribute
	47	!
[1323]	48	! 1322 2014-03-20 16:38:49Z raasch
	49	! REAL constants defined as wp-kind
	50	!
[1321]	51	! 1320 2014-03-20 08:40:49Z raasch
[1320]	52	! ONLY-attribute added to USE-statements,
	53	! kind-parameters added to all INTEGER and REAL declaration statements,
	54	! kinds are defined in new module kinds,
	55	! comment fields (!:) to be used for variable explanations added to
	56	! all variable declaration statements
[1000]	57	!
[1242]	58	! 1241 2013-10-30 11:36:58Z heinze
	59	! hyp and rho have to be calculated at each time step if data from external
	60	! file LSF_DATA are used
	61	!
[1116]	62	! 1115 2013-03-26 18:16:16Z hoffmann
	63	! microphyical tendencies are calculated in microphysics_control in an optimized
	64	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
	65	!
[1107]	66	! 1106 2013-03-04 05:31:38Z raasch
	67	! small changes in code formatting
	68	!
[1093]	69	! 1092 2013-02-02 11:24:22Z raasch
	70	! unused variables removed
	71	! file put under GPL
	72	!
[1066]	73	! 1065 2012-11-22 17:42:36Z hoffmann
	74	! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]	75	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]	76	! and Stevens, 2010).
	77	!
[1054]	78	! 1053 2012-11-13 17:11:03Z hoffmann
	79	! initial revision
[1000]	80	!
	81	! Description:
	82	! ------------
[1682]	83	!> Calculate cloud microphysics according to the two moment bulk
	84	!> scheme by Seifert and Beheng (2006).
[1000]	85	!------------------------------------------------------------------------------!
[1682]	86	MODULE microphysics_mod
	87
[1000]	88
	89	PRIVATE
[1115]	90	PUBLIC microphysics_control
[1000]	91
[1115]	92	INTERFACE microphysics_control
	93	MODULE PROCEDURE microphysics_control
	94	MODULE PROCEDURE microphysics_control_ij
	95	END INTERFACE microphysics_control
[1022]	96
[1115]	97	INTERFACE adjust_cloud
	98	MODULE PROCEDURE adjust_cloud
	99	MODULE PROCEDURE adjust_cloud_ij
	100	END INTERFACE adjust_cloud
	101
[1000]	102	INTERFACE autoconversion
	103	MODULE PROCEDURE autoconversion
	104	MODULE PROCEDURE autoconversion_ij
	105	END INTERFACE autoconversion
	106
	107	INTERFACE accretion
	108	MODULE PROCEDURE accretion
	109	MODULE PROCEDURE accretion_ij
	110	END INTERFACE accretion
[1005]	111
	112	INTERFACE selfcollection_breakup
	113	MODULE PROCEDURE selfcollection_breakup
	114	MODULE PROCEDURE selfcollection_breakup_ij
	115	END INTERFACE selfcollection_breakup
[1012]	116
	117	INTERFACE evaporation_rain
	118	MODULE PROCEDURE evaporation_rain
	119	MODULE PROCEDURE evaporation_rain_ij
	120	END INTERFACE evaporation_rain
	121
	122	INTERFACE sedimentation_cloud
	123	MODULE PROCEDURE sedimentation_cloud
	124	MODULE PROCEDURE sedimentation_cloud_ij
	125	END INTERFACE sedimentation_cloud
[1000]	126
[1012]	127	INTERFACE sedimentation_rain
	128	MODULE PROCEDURE sedimentation_rain
	129	MODULE PROCEDURE sedimentation_rain_ij
	130	END INTERFACE sedimentation_rain
	131
[1000]	132	CONTAINS
	133
	134
	135	!------------------------------------------------------------------------------!
[1682]	136	! Description:
	137	! ------------
	138	!> Call for all grid points
[1000]	139	!------------------------------------------------------------------------------!
[1115]	140	SUBROUTINE microphysics_control
[1022]	141
[1361]	142	USE arrays_3d, &
	143	ONLY: hyp, nr, pt, pt_init, q, qc, qr, zu
	144
	145	USE cloud_parameters, &
	146	ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
	147
	148	USE control_parameters, &
	149	ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
	150	g, intermediate_timestep_count, &
	151	large_scale_forcing, lsf_surf, precipitation, pt_surface, &
	152	rho_surface,surface_pressure
	153
	154	USE indices, &
	155	ONLY: nzb, nzt
	156
[1320]	157	USE kinds
[1115]	158
[1361]	159	USE statistics, &
	160	ONLY: weight_pres
	161
[1115]	162	IMPLICIT NONE
	163
[1682]	164	INTEGER(iwp) :: i !<
	165	INTEGER(iwp) :: j !<
	166	INTEGER(iwp) :: k !<
[1115]	167
[1682]	168	REAL(wp) :: t_surface !<
[1361]	169
	170	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	171	!
	172	!-- Calculate:
	173	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	174	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	175	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
	176	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
	177	DO k = nzb, nzt+1
	178	hyp(k) = surface_pressure * 100.0_wp * &
	179	( ( t_surface - g / cp * zu(k) ) / &
	180	t_surface )**(1.0_wp / 0.286_wp)
	181	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	182	t_d_pt(k) = 1.0_wp / pt_d_t(k)
	183	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1115]	184	ENDDO
[1361]	185	!
	186	!-- Compute reference density
	187	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
	188	ENDIF
[1115]	189
[1361]	190	!
	191	!-- Compute length of time step
	192	IF ( call_microphysics_at_all_substeps ) THEN
	193	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	194	ELSE
	195	dt_micro = dt_3d
	196	ENDIF
	197
	198	!
	199	!-- Compute cloud physics
	200	IF ( precipitation ) THEN
	201	CALL adjust_cloud
	202	CALL autoconversion
	203	CALL accretion
	204	CALL selfcollection_breakup
	205	CALL evaporation_rain
	206	CALL sedimentation_rain
	207	ENDIF
	208
	209	IF ( drizzle ) CALL sedimentation_cloud
	210
[1115]	211	END SUBROUTINE microphysics_control
	212
[1682]	213	!------------------------------------------------------------------------------!
	214	! Description:
	215	! ------------
	216	!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
	217	!> evaporation of rain drops due to too small or too big weights
	218	!> of rain drops (Stevens and Seifert, 2008).
	219	!------------------------------------------------------------------------------!
[1115]	220	SUBROUTINE adjust_cloud
	221
[1361]	222	USE arrays_3d, &
	223	ONLY: qr, nr
	224
	225	USE cloud_parameters, &
	226	ONLY: eps_sb, xrmin, xrmax, hyrho, k_cc, x0
	227
	228	USE cpulog, &
	229	ONLY: cpu_log, log_point_s
	230
	231	USE indices, &
	232	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	233
[1320]	234	USE kinds
[1022]	235
	236	IMPLICIT NONE
	237
[1682]	238	INTEGER(iwp) :: i !<
	239	INTEGER(iwp) :: j !<
	240	INTEGER(iwp) :: k !<
[1022]	241
[1361]	242	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
	243
[1022]	244	DO i = nxl, nxr
	245	DO j = nys, nyn
[1115]	246	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	247	IF ( qr(k,j,i) <= eps_sb ) THEN
	248	qr(k,j,i) = 0.0_wp
	249	nr(k,j,i) = 0.0_wp
	250	ELSE
	251	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
	252	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
	253	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
	254	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
	255	ENDIF
	256	ENDIF
[1022]	257	ENDDO
	258	ENDDO
	259	ENDDO
	260
[1361]	261	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
	262
[1115]	263	END SUBROUTINE adjust_cloud
[1022]	264
[1106]	265
[1682]	266	!------------------------------------------------------------------------------!
	267	! Description:
	268	! ------------
	269	!> Autoconversion rate (Seifert and Beheng, 2006).
	270	!------------------------------------------------------------------------------!
[1000]	271	SUBROUTINE autoconversion
	272
[1361]	273	USE arrays_3d, &
	274	ONLY: diss, dzu, nr, qc, qr
	275
	276	USE cloud_parameters, &
	277	ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, &
	278	c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, &
	279	nc_const, x0
	280
	281	USE control_parameters, &
	282	ONLY: dt_micro, rho_surface, turbulence
	283
	284	USE cpulog, &
	285	ONLY: cpu_log, log_point_s
	286
	287	USE grid_variables, &
	288	ONLY: dx, dy
	289
	290	USE indices, &
	291	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	292
[1320]	293	USE kinds
[1000]	294
	295	IMPLICIT NONE
	296
[1682]	297	INTEGER(iwp) :: i !<
	298	INTEGER(iwp) :: j !<
	299	INTEGER(iwp) :: k !<
[1000]	300
[1682]	301	REAL(wp) :: alpha_cc !<
	302	REAL(wp) :: autocon !<
	303	REAL(wp) :: dissipation !<
	304	REAL(wp) :: k_au !<
	305	REAL(wp) :: l_mix !<
	306	REAL(wp) :: nu_c !<
	307	REAL(wp) :: phi_au !<
	308	REAL(wp) :: r_cc !<
	309	REAL(wp) :: rc !<
	310	REAL(wp) :: re_lambda !<
	311	REAL(wp) :: selfcoll !<
	312	REAL(wp) :: sigma_cc !<
	313	REAL(wp) :: tau_cloud !<
	314	REAL(wp) :: xc !<
[1361]	315
	316	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
	317
[1000]	318	DO i = nxl, nxr
	319	DO j = nys, nyn
[1115]	320	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	321
[1361]	322	IF ( qc(k,j,i) > eps_sb ) THEN
	323
	324	k_au = k_cc / ( 20.0_wp * x0 )
	325	!
	326	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	327	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
	328	tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
	329	!
	330	!-- Universal function for autoconversion process
	331	!-- (Seifert and Beheng, 2006):
	332	phi_au = 600.0_wp * tau_cloud*0.68_wp &
	333	( 1.0_wp - tau_cloud0.68_wp )3
	334	!
	335	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
	336	!-- (Use constant nu_c = 1.0_wp instead?)
	337	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
	338	!
	339	!-- Mean weight of cloud droplets:
	340	xc = hyrho(k) * qc(k,j,i) / nc_const
	341	!
	342	!-- Parameterized turbulence effects on autoconversion (Seifert,
	343	!-- Nuijens and Stevens, 2010)
	344	IF ( turbulence ) THEN
	345	!
	346	!-- Weight averaged radius of cloud droplets:
	347	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	348
	349	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	350	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	351	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
	352	!
	353	!-- Mixing length (neglecting distance to ground and
	354	!-- stratification)
	355	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
	356	!
	357	!-- Limit dissipation rate according to Seifert, Nuijens and
	358	!-- Stevens (2010)
	359	dissipation = MIN( 0.06_wp, diss(k,j,i) )
	360	!
	361	!-- Compute Taylor-microscale Reynolds number:
	362	re_lambda = 6.0_wp / 11.0_wp * &
	363	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	364	SQRT( 15.0_wp / kin_vis_air ) * &
	365	dissipation**( 1.0_wp / 6.0_wp )
	366	!
	367	!-- The factor of 1.0E4 is needed to convert the dissipation
	368	!-- rate from m2 s-3 to cm2 s-3.
	369	k_au = k_au * ( 1.0_wp + &
	370	dissipation * 1.0E4_wp * &
	371	( re_lambda * 1.0E-3_wp )*0.25_wp &
	372	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
	373	r_cc ) / &
	374	sigma_cc )**2 &
	375	) + beta_cc &
	376	) &
	377	)
	378	ENDIF
	379	!
	380	!-- Autoconversion rate (Seifert and Beheng, 2006):
	381	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	382	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
	383	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
	384	rho_surface
	385	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
	386
	387	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
	388	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
	389	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
	390
	391	ENDIF
	392
[1000]	393	ENDDO
	394	ENDDO
	395	ENDDO
	396
[1361]	397	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
	398
[1000]	399	END SUBROUTINE autoconversion
	400
[1106]	401
[1682]	402	!------------------------------------------------------------------------------!
	403	! Description:
	404	! ------------
	405	!> Accretion rate (Seifert and Beheng, 2006).
	406	!------------------------------------------------------------------------------!
[1005]	407	SUBROUTINE accretion
[1000]	408
[1361]	409	USE arrays_3d, &
	410	ONLY: diss, qc, qr
	411
	412	USE cloud_parameters, &
	413	ONLY: eps_sb, hyrho, k_cr0
	414
	415	USE control_parameters, &
	416	ONLY: dt_micro, rho_surface, turbulence
	417
	418	USE cpulog, &
	419	ONLY: cpu_log, log_point_s
	420
	421	USE indices, &
	422	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	423
[1320]	424	USE kinds
[1005]	425
[1000]	426	IMPLICIT NONE
	427
[1682]	428	INTEGER(iwp) :: i !<
	429	INTEGER(iwp) :: j !<
	430	INTEGER(iwp) :: k !<
[1000]	431
[1682]	432	REAL(wp) :: accr !<
	433	REAL(wp) :: k_cr !<
	434	REAL(wp) :: phi_ac !<
	435	REAL(wp) :: tau_cloud !<
	436	REAL(wp) :: xc !<
[1361]	437
	438	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
	439
[1005]	440	DO i = nxl, nxr
	441	DO j = nys, nyn
[1115]	442	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	443
[1361]	444	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
	445	!
	446	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	447	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
	448	!
	449	!-- Universal function for accretion process (Seifert and
	450	!-- Beheng, 2001):
	451	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
	452	!
	453	!-- Parameterized turbulence effects on autoconversion (Seifert,
	454	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
	455	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
	456	IF ( turbulence ) THEN
	457	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	458	MIN( 600.0_wp, &
	459	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	460	)
	461	ELSE
	462	k_cr = k_cr0
	463	ENDIF
	464	!
	465	!-- Accretion rate (Seifert and Beheng, 2006):
	466	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
	467	SQRT( rho_surface * hyrho(k) )
	468	accr = MIN( accr, qc(k,j,i) / dt_micro )
	469
	470	qr(k,j,i) = qr(k,j,i) + accr * dt_micro
	471	qc(k,j,i) = qc(k,j,i) - accr * dt_micro
	472
	473	ENDIF
	474
[1005]	475	ENDDO
	476	ENDDO
[1000]	477	ENDDO
	478
[1361]	479	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
	480
[1005]	481	END SUBROUTINE accretion
[1000]	482
[1106]	483
[1682]	484	!------------------------------------------------------------------------------!
	485	! Description:
	486	! ------------
	487	!> Collisional breakup rate (Seifert, 2008).
	488	!------------------------------------------------------------------------------!
[1005]	489	SUBROUTINE selfcollection_breakup
[1000]	490
[1361]	491	USE arrays_3d, &
	492	ONLY: nr, qr
	493
	494	USE cloud_parameters, &
	495	ONLY: dpirho_l, eps_sb, hyrho, k_br, k_rr
	496
	497	USE control_parameters, &
	498	ONLY: dt_micro, rho_surface
	499
	500	USE cpulog, &
	501	ONLY: cpu_log, log_point_s
	502
	503	USE indices, &
	504	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	505
[1320]	506	USE kinds
[1361]	507
[1000]	508	IMPLICIT NONE
	509
[1682]	510	INTEGER(iwp) :: i !<
	511	INTEGER(iwp) :: j !<
	512	INTEGER(iwp) :: k !<
[1000]	513
[1682]	514	REAL(wp) :: breakup !<
	515	REAL(wp) :: dr !<
	516	REAL(wp) :: phi_br !<
	517	REAL(wp) :: selfcoll !<
[1361]	518
	519	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
	520
[1000]	521	DO i = nxl, nxr
	522	DO j = nys, nyn
[1115]	523	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	524	IF ( qr(k,j,i) > eps_sb ) THEN
	525	!
	526	!-- Selfcollection rate (Seifert and Beheng, 2001):
	527	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
	528	SQRT( hyrho(k) * rho_surface )
	529	!
	530	!-- Weight averaged diameter of rain drops:
	531	dr = ( hyrho(k) * qr(k,j,i) / &
	532	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	533	!
	534	!-- Collisional breakup rate (Seifert, 2008):
	535	IF ( dr >= 0.3E-3_wp ) THEN
	536	phi_br = k_br * ( dr - 1.1E-3_wp )
	537	breakup = selfcoll * ( phi_br + 1.0_wp )
	538	ELSE
	539	breakup = 0.0_wp
	540	ENDIF
[1000]	541
[1361]	542	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
	543	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
	544
	545	ENDIF
[1000]	546	ENDDO
	547	ENDDO
	548	ENDDO
	549
[1361]	550	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
	551
[1005]	552	END SUBROUTINE selfcollection_breakup
[1000]	553
[1106]	554
[1682]	555	!------------------------------------------------------------------------------!
	556	! Description:
	557	! ------------
	558	!> Evaporation of precipitable water. Condensation is neglected for
	559	!> precipitable water.
	560	!------------------------------------------------------------------------------!
[1012]	561	SUBROUTINE evaporation_rain
[1000]	562
[1361]	563	USE arrays_3d, &
	564	ONLY: hyp, nr, pt, q, qc, qr
	565
	566	USE cloud_parameters, &
	567	ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
	568	dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r, &
	569	l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l, &
	570	t_d_pt, ventilation_effect
	571
	572	USE constants, &
	573	ONLY: pi
	574
	575	USE control_parameters, &
	576	ONLY: dt_micro
	577
	578	USE cpulog, &
	579	ONLY: cpu_log, log_point_s
	580
	581	USE indices, &
	582	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	583
[1320]	584	USE kinds
[1012]	585
	586	IMPLICIT NONE
	587
[1682]	588	INTEGER(iwp) :: i !<
	589	INTEGER(iwp) :: j !<
	590	INTEGER(iwp) :: k !<
[1361]	591
[1682]	592	REAL(wp) :: alpha !<
	593	REAL(wp) :: dr !<
	594	REAL(wp) :: e_s !<
	595	REAL(wp) :: evap !<
	596	REAL(wp) :: evap_nr !<
	597	REAL(wp) :: f_vent !<
	598	REAL(wp) :: g_evap !<
	599	REAL(wp) :: lambda_r !<
	600	REAL(wp) :: mu_r !<
	601	REAL(wp) :: mu_r_2 !<
	602	REAL(wp) :: mu_r_5d2 !<
	603	REAL(wp) :: nr_0 !<
	604	REAL(wp) :: q_s !<
	605	REAL(wp) :: sat !<
	606	REAL(wp) :: t_l !<
	607	REAL(wp) :: temp !<
	608	REAL(wp) :: xr !<
[1361]	609
	610	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
	611
[1012]	612	DO i = nxl, nxr
	613	DO j = nys, nyn
[1115]	614	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	615	IF ( qr(k,j,i) > eps_sb ) THEN
	616	!
	617	!-- Actual liquid water temperature:
	618	t_l = t_d_pt(k) * pt(k,j,i)
	619	!
	620	!-- Saturation vapor pressure at t_l:
	621	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	622	( t_l - 35.86_wp ) &
	623	)
	624	!
	625	!-- Computation of saturation humidity:
	626	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	627	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	628	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
	629	( 1.0_wp + alpha * q_s )
	630	!
	631	!-- Supersaturation:
	632	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
	633	!
	634	!-- Evaporation needs only to be calculated in subsaturated regions
	635	IF ( sat < 0.0_wp ) THEN
	636	!
	637	!-- Actual temperature:
	638	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
	639
	640	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	641	l_v / ( thermal_conductivity_l * temp ) &
	642	+ r_v * temp / ( diff_coeff_l * e_s ) &
	643	)
	644	!
	645	!-- Mean weight of rain drops
	646	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
	647	!
	648	!-- Weight averaged diameter of rain drops:
	649	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
	650	!
	651	!-- Compute ventilation factor and intercept parameter
	652	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	653	IF ( ventilation_effect ) THEN
	654	!
	655	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
	656	!-- 2005; Stevens and Seifert, 2008):
	657	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	658	( dr - 1.4E-3_wp ) ) )
	659	!
	660	!-- Slope parameter of gamma distribution (Seifert, 2008):
	661	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	662	( mu_r + 1.0_wp ) &
	663	)**( 1.0_wp / 3.0_wp ) / dr
[1012]	664
[1361]	665	mu_r_2 = mu_r + 2.0_wp
	666	mu_r_5d2 = mu_r + 2.5_wp
	667
	668	f_vent = a_vent * gamm( mu_r_2 ) * &
	669	lambda_r*( -mu_r_2 ) + b_vent &
	670	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
	671	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
	672	( 1.0_wp - &
	673	0.5_wp * ( b_term / a_term ) * &
	674	( lambda_r / ( c_term + lambda_r ) &
	675	)**mu_r_5d2 - &
	676	0.125_wp * ( b_term / a_term )*2 &
	677	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	678	)**mu_r_5d2 - &
	679	0.0625_wp * ( b_term / a_term )*3 &
	680	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	681	)**mu_r_5d2 - &
	682	0.0390625_wp * ( b_term / a_term )*4 &
	683	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	684	)**mu_r_5d2 &
	685	)
	686
	687	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
	688	gamm( mu_r + 1.0_wp )
	689	ELSE
	690	f_vent = 1.0_wp
	691	nr_0 = nr(k,j,i) * dr
	692	ENDIF
	693	!
	694	!-- Evaporation rate of rain water content (Seifert and
	695	!-- Beheng, 2006):
	696	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
	697	hyrho(k)
	698	evap = MAX( evap, -qr(k,j,i) / dt_micro )
	699	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	700	-nr(k,j,i) / dt_micro )
	701
	702	qr(k,j,i) = qr(k,j,i) + evap * dt_micro
	703	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
	704
	705	ENDIF
	706	ENDIF
	707
[1012]	708	ENDDO
	709	ENDDO
	710	ENDDO
	711
[1361]	712	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
	713
[1012]	714	END SUBROUTINE evaporation_rain
	715
[1106]	716
[1682]	717	!------------------------------------------------------------------------------!
	718	! Description:
	719	! ------------
	720	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
	721	!------------------------------------------------------------------------------!
[1012]	722	SUBROUTINE sedimentation_cloud
	723
[1361]	724	USE arrays_3d, &
	725	ONLY: ddzu, dzu, pt, q, qc
	726
	727	USE cloud_parameters, &
	728	ONLY: eps_sb, hyrho, l_d_cp, nc_const, pt_d_t, sed_qc_const
	729
	730	USE constants, &
	731	ONLY: pi
	732
	733	USE control_parameters, &
	734	ONLY: dt_do2d_xy, dt_micro, intermediate_timestep_count
	735
	736	USE cpulog, &
	737	ONLY: cpu_log, log_point_s
	738
	739	USE indices, &
	740	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	741
[1320]	742	USE kinds
[1361]	743
[1012]	744	IMPLICIT NONE
	745
[1682]	746	INTEGER(iwp) :: i !<
	747	INTEGER(iwp) :: j !<
	748	INTEGER(iwp) :: k !<
[1361]	749
[1682]	750	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1361]	751
	752	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
	753
	754	sed_qc(nzt+1) = 0.0_wp
	755
[1012]	756	DO i = nxl, nxr
	757	DO j = nys, nyn
[1361]	758	DO k = nzt, nzb_s_inner(j,i)+1, -1
[1012]	759
[1361]	760	IF ( qc(k,j,i) > eps_sb ) THEN
	761	sed_qc(k) = sed_qc_const * nc_const*( -2.0_wp / 3.0_wp ) &
	762	( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
	763	ELSE
	764	sed_qc(k) = 0.0_wp
	765	ENDIF
	766
	767	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
	768	dt_micro + sed_qc(k+1) &
	769	)
	770
	771	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
	772	ddzu(k+1) / hyrho(k) * dt_micro
	773	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
	774	ddzu(k+1) / hyrho(k) * dt_micro
	775	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
	776	ddzu(k+1) / hyrho(k) * l_d_cp * &
	777	pt_d_t(k) * dt_micro
	778
[1012]	779	ENDDO
	780	ENDDO
	781	ENDDO
	782
[1361]	783	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
	784
[1012]	785	END SUBROUTINE sedimentation_cloud
	786
[1106]	787
[1682]	788	!------------------------------------------------------------------------------!
	789	! Description:
	790	! ------------
	791	!> Computation of sedimentation flux. Implementation according to Stevens
	792	!> and Seifert (2008). Code is based on UCLA-LES.
	793	!------------------------------------------------------------------------------!
[1012]	794	SUBROUTINE sedimentation_rain
	795
[1361]	796	USE arrays_3d, &
	797	ONLY: ddzu, dzu, nr, pt, q, qr
	798
	799	USE cloud_parameters, &
	800	ONLY: a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho, &
	801	limiter_sedimentation, l_d_cp, precipitation_amount, prr, &
	802	pt_d_t, stp
	803
	804	USE control_parameters, &
	805	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro, &
	806	dt_3d, intermediate_timestep_count, &
	807	intermediate_timestep_count_max, &
	808	precipitation_amount_interval, time_do2d_xy
	809
	810	USE cpulog, &
	811	ONLY: cpu_log, log_point_s
	812
	813	USE indices, &
	814	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	815
[1320]	816	USE kinds
[1012]	817
[1361]	818	USE statistics, &
	819	ONLY: weight_substep
	820
[1012]	821	IMPLICIT NONE
	822
[1682]	823	INTEGER(iwp) :: i !<
	824	INTEGER(iwp) :: j !<
	825	INTEGER(iwp) :: k !<
	826	INTEGER(iwp) :: k_run !<
[1361]	827
[1682]	828	REAL(wp) :: c_run !<
	829	REAL(wp) :: d_max !<
	830	REAL(wp) :: d_mean !<
	831	REAL(wp) :: d_min !<
	832	REAL(wp) :: dr !<
	833	REAL(wp) :: dt_sedi !<
	834	REAL(wp) :: flux !<
	835	REAL(wp) :: lambda_r !<
	836	REAL(wp) :: mu_r !<
	837	REAL(wp) :: z_run !<
[1361]	838
[1682]	839	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	840	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	841	REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !<
	842	REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !<
	843	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	844	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	845	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	846	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	847	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	848	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1361]	849
	850	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
[1682]	851
[1361]	852	IF ( intermediate_timestep_count == 1 ) prr(:,:,:) = 0.0_wp
	853	!
	854	!-- Compute velocities
[1012]	855	DO i = nxl, nxr
	856	DO j = nys, nyn
[1115]	857	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	858	IF ( qr(k,j,i) > eps_sb ) THEN
	859	!
	860	!-- Weight averaged diameter of rain drops:
	861	dr = ( hyrho(k) * qr(k,j,i) / &
	862	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	863	!
	864	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	865	!-- Stevens and Seifert, 2008):
	866	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	867	( dr - 1.4E-3_wp ) ) )
	868	!
	869	!-- Slope parameter of gamma distribution (Seifert, 2008):
	870	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	871	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1012]	872
[1361]	873	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	874	a_term - b_term * ( 1.0_wp + &
	875	c_term / &
	876	lambda_r )*( -1.0_wp &
	877	( mu_r + 1.0_wp ) ) &
	878	) &
	879	)
	880
	881	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	882	a_term - b_term * ( 1.0_wp + &
	883	c_term / &
	884	lambda_r )*( -1.0_wp &
	885	( mu_r + 4.0_wp ) ) &
	886	) &
	887	)
	888	ELSE
	889	w_nr(k) = 0.0_wp
	890	w_qr(k) = 0.0_wp
	891	ENDIF
[1012]	892	ENDDO
[1361]	893	!
	894	!-- Adjust boundary values
	895	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
	896	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
	897	w_nr(nzt+1) = 0.0_wp
	898	w_qr(nzt+1) = 0.0_wp
	899	!
	900	!-- Compute Courant number
	901	DO k = nzb_s_inner(j,i)+1, nzt
	902	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
	903	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
	904	dt_micro * ddzu(k)
	905	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
	906	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
	907	dt_micro * ddzu(k)
	908	ENDDO
	909	!
	910	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	911	IF ( limiter_sedimentation ) THEN
	912
	913	DO k = nzb_s_inner(j,i)+1, nzt
[1646]	914	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
[1361]	915	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
	916	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
	917
	918	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	919	2.0_wp * d_max, &
	920	ABS( d_mean ) )
	921
[1646]	922	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
[1361]	923	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
	924	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
	925
	926	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	927	2.0_wp * d_max, &
	928	ABS( d_mean ) )
	929	ENDDO
	930
	931	ELSE
	932
	933	nr_slope = 0.0_wp
	934	qr_slope = 0.0_wp
	935
	936	ENDIF
	937
	938	sed_nr(nzt+1) = 0.0_wp
	939	sed_qr(nzt+1) = 0.0_wp
	940	!
	941	!-- Compute sedimentation flux
	942	DO k = nzt, nzb_s_inner(j,i)+1, -1
	943	!
	944	!-- Sum up all rain drop number densities which contribute to the flux
	945	!-- through k-1/2
	946	flux = 0.0_wp
	947	z_run = 0.0_wp ! height above z(k)
	948	k_run = k
	949	c_run = MIN( 1.0_wp, c_nr(k) )
	950	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	951	flux = flux + hyrho(k_run) * &
	952	( nr(k_run,j,i) + nr_slope(k_run) * &
	953	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
	954	z_run = z_run + dzu(k_run)
	955	k_run = k_run + 1
	956	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
	957	ENDDO
	958	!
	959	!-- It is not allowed to sediment more rain drop number density than
	960	!-- available
	961	flux = MIN( flux, &
	962	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
	963	dt_micro &
	964	)
	965
	966	sed_nr(k) = flux / dt_micro
	967	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
	968	ddzu(k+1) / hyrho(k) * dt_micro
	969	!
	970	!-- Sum up all rain water content which contributes to the flux
	971	!-- through k-1/2
	972	flux = 0.0_wp
	973	z_run = 0.0_wp ! height above z(k)
	974	k_run = k
	975	c_run = MIN( 1.0_wp, c_qr(k) )
	976
	977	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	978
	979	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
	980	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
	981	0.5_wp ) * c_run * dzu(k_run)
	982	z_run = z_run + dzu(k_run)
	983	k_run = k_run + 1
	984	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
	985
	986	ENDDO
	987	!
	988	!-- It is not allowed to sediment more rain water content than
	989	!-- available
	990	flux = MIN( flux, &
	991	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
	992	dt_micro &
	993	)
	994
	995	sed_qr(k) = flux / dt_micro
	996
	997	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
	998	ddzu(k+1) / hyrho(k) * dt_micro
	999	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
	1000	ddzu(k+1) / hyrho(k) * dt_micro
	1001	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
	1002	ddzu(k+1) / hyrho(k) * l_d_cp * &
	1003	pt_d_t(k) * dt_micro
	1004	!
	1005	!-- Compute the rain rate
	1006	IF ( call_microphysics_at_all_substeps ) THEN
	1007	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
	1008	weight_substep(intermediate_timestep_count)
	1009	ELSE
	1010	prr(k,j,i) = sed_qr(k) / hyrho(k)
	1011	ENDIF
	1012
	1013	ENDDO
[1012]	1014	ENDDO
	1015	ENDDO
	1016
[1361]	1017	!
	1018	!-- Precipitation amount
	1019	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
	1020	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
	1021	precipitation_amount_interval ) THEN
	1022	DO i = nxl, nxr
	1023	DO j = nys, nyn
	1024	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	1025	prr(nzb_s_inner(j,i)+1,j,i) * &
	1026	hyrho(nzb_s_inner(j,i)+1) * dt_3d
	1027	ENDDO
	1028	ENDDO
	1029	ENDIF
	1030
	1031	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
	1032
[1012]	1033	END SUBROUTINE sedimentation_rain
	1034
	1035
[1000]	1036	!------------------------------------------------------------------------------!
[1682]	1037	! Description:
	1038	! ------------
	1039	!> Call for grid point i,j
[1000]	1040	!------------------------------------------------------------------------------!
[1022]	1041
[1115]	1042	SUBROUTINE microphysics_control_ij( i, j )
	1043
[1320]	1044	USE arrays_3d, &
[1361]	1045	ONLY: hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc, &
	1046	qc_1d, qr, qr_1d, zu
[1115]	1047
[1320]	1048	USE cloud_parameters, &
	1049	ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
	1050
	1051	USE control_parameters, &
[1361]	1052	ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
	1053	g, intermediate_timestep_count, large_scale_forcing, &
	1054	lsf_surf, precipitation, pt_surface, &
[1320]	1055	rho_surface,surface_pressure
	1056
	1057	USE indices, &
	1058	ONLY: nzb, nzt
	1059
	1060	USE kinds
	1061
	1062	USE statistics, &
	1063	ONLY: weight_pres
	1064
[1022]	1065	IMPLICIT NONE
	1066
[1682]	1067	INTEGER(iwp) :: i !<
	1068	INTEGER(iwp) :: j !<
	1069	INTEGER(iwp) :: k !<
[1115]	1070
[1682]	1071	REAL(wp) :: t_surface !<
[1320]	1072
[1361]	1073	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[1241]	1074	!
	1075	!-- Calculate:
	1076	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	1077	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	1078	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1353]	1079	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]	1080	DO k = nzb, nzt+1
[1353]	1081	hyp(k) = surface_pressure * 100.0_wp * &
[1361]	1082	( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
[1353]	1083	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	1084	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[1241]	1085	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
	1086	ENDDO
	1087	!
	1088	!-- Compute reference density
[1353]	1089	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]	1090	ENDIF
	1091
[1361]	1092	!
	1093	!-- Compute length of time step
	1094	IF ( call_microphysics_at_all_substeps ) THEN
	1095	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	1096	ELSE
	1097	dt_micro = dt_3d
	1098	ENDIF
[1241]	1099
[1115]	1100	!
[1361]	1101	!-- Use 1d arrays
[1115]	1102	q_1d(:) = q(:,j,i)
	1103	pt_1d(:) = pt(:,j,i)
	1104	qc_1d(:) = qc(:,j,i)
	1105	nc_1d(:) = nc_const
	1106	IF ( precipitation ) THEN
	1107	qr_1d(:) = qr(:,j,i)
	1108	nr_1d(:) = nr(:,j,i)
	1109	ENDIF
[1361]	1110
[1115]	1111	!
	1112	!-- Compute cloud physics
	1113	IF ( precipitation ) THEN
[1361]	1114	CALL adjust_cloud( i,j )
[1115]	1115	CALL autoconversion( i,j )
	1116	CALL accretion( i,j )
	1117	CALL selfcollection_breakup( i,j )
	1118	CALL evaporation_rain( i,j )
	1119	CALL sedimentation_rain( i,j )
	1120	ENDIF
	1121
	1122	IF ( drizzle ) CALL sedimentation_cloud( i,j )
[1361]	1123
[1115]	1124	!
[1361]	1125	!-- Store results on the 3d arrays
	1126	q(:,j,i) = q_1d(:)
	1127	pt(:,j,i) = pt_1d(:)
[1115]	1128	IF ( precipitation ) THEN
[1361]	1129	qr(:,j,i) = qr_1d(:)
	1130	nr(:,j,i) = nr_1d(:)
[1115]	1131	ENDIF
	1132
	1133	END SUBROUTINE microphysics_control_ij
	1134
[1682]	1135	!------------------------------------------------------------------------------!
	1136	! Description:
	1137	! ------------
	1138	!> Adjust number of raindrops to avoid nonlinear effects in
	1139	!> sedimentation and evaporation of rain drops due to too small or
	1140	!> too big weights of rain drops (Stevens and Seifert, 2008).
	1141	!> The same procedure is applied to cloud droplets if they are determined
	1142	!> prognostically. Call for grid point i,j
	1143	!------------------------------------------------------------------------------!
[1115]	1144	SUBROUTINE adjust_cloud_ij( i, j )
	1145
[1320]	1146	USE arrays_3d, &
[1361]	1147	ONLY: qr_1d, nr_1d
[1115]	1148
[1320]	1149	USE cloud_parameters, &
	1150	ONLY: eps_sb, xrmin, xrmax, hyrho, k_cc, x0
	1151
	1152	USE indices, &
	1153	ONLY: nzb, nzb_s_inner, nzt
	1154
	1155	USE kinds
	1156
[1115]	1157	IMPLICIT NONE
	1158
[1682]	1159	INTEGER(iwp) :: i !<
	1160	INTEGER(iwp) :: j !<
	1161	INTEGER(iwp) :: k !<
	1162
[1115]	1163	DO k = nzb_s_inner(j,i)+1, nzt
[1022]	1164
[1361]	1165	IF ( qr_1d(k) <= eps_sb ) THEN
	1166	qr_1d(k) = 0.0_wp
	1167	nr_1d(k) = 0.0_wp
[1065]	1168	ELSE
[1022]	1169	!
[1048]	1170	!-- Adjust number of raindrops to avoid nonlinear effects in
	1171	!-- sedimentation and evaporation of rain drops due to too small or
[1065]	1172	!-- too big weights of rain drops (Stevens and Seifert, 2008).
[1361]	1173	IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN
	1174	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
	1175	ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN
	1176	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
[1048]	1177	ENDIF
[1115]	1178
[1022]	1179	ENDIF
[1115]	1180
[1022]	1181	ENDDO
	1182
[1115]	1183	END SUBROUTINE adjust_cloud_ij
[1022]	1184
[1106]	1185
[1682]	1186	!------------------------------------------------------------------------------!
	1187	! Description:
	1188	! ------------
	1189	!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
	1190	!------------------------------------------------------------------------------!
[1005]	1191	SUBROUTINE autoconversion_ij( i, j )
[1000]	1192
[1320]	1193	USE arrays_3d, &
	1194	ONLY: diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
[1115]	1195
[1320]	1196	USE cloud_parameters, &
	1197	ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, &
	1198	c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
	1199
	1200	USE control_parameters, &
	1201	ONLY: dt_micro, rho_surface, turbulence
	1202
	1203	USE grid_variables, &
	1204	ONLY: dx, dy
	1205
	1206	USE indices, &
	1207	ONLY: nzb, nzb_s_inner, nzt
	1208
	1209	USE kinds
	1210
[1000]	1211	IMPLICIT NONE
	1212
[1682]	1213	INTEGER(iwp) :: i !<
	1214	INTEGER(iwp) :: j !<
	1215	INTEGER(iwp) :: k !<
[1000]	1216
[1682]	1217	REAL(wp) :: alpha_cc !<
	1218	REAL(wp) :: autocon !<
	1219	REAL(wp) :: dissipation !<
	1220	REAL(wp) :: k_au !<
	1221	REAL(wp) :: l_mix !<
	1222	REAL(wp) :: nu_c !<
	1223	REAL(wp) :: phi_au !<
	1224	REAL(wp) :: r_cc !<
	1225	REAL(wp) :: rc !<
	1226	REAL(wp) :: re_lambda !<
	1227	REAL(wp) :: selfcoll !<
	1228	REAL(wp) :: sigma_cc !<
	1229	REAL(wp) :: tau_cloud !<
	1230	REAL(wp) :: xc !<
[1106]	1231
[1115]	1232	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	1233
[1115]	1234	IF ( qc_1d(k) > eps_sb ) THEN
[1361]	1235
	1236	k_au = k_cc / ( 20.0_wp * x0 )
[1012]	1237	!
[1048]	1238	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	1239	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
	1240	tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
[1012]	1241	!
	1242	!-- Universal function for autoconversion process
	1243	!-- (Seifert and Beheng, 2006):
[1361]	1244	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
[1012]	1245	!
	1246	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]	1247	!-- (Use constant nu_c = 1.0_wp instead?)
[1361]	1248	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
[1012]	1249	!
	1250	!-- Mean weight of cloud droplets:
[1115]	1251	xc = hyrho(k) * qc_1d(k) / nc_1d(k)
[1012]	1252	!
[1065]	1253	!-- Parameterized turbulence effects on autoconversion (Seifert,
	1254	!-- Nuijens and Stevens, 2010)
	1255	IF ( turbulence ) THEN
	1256	!
	1257	!-- Weight averaged radius of cloud droplets:
[1353]	1258	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]	1259
[1353]	1260	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	1261	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	1262	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
[1065]	1263	!
	1264	!-- Mixing length (neglecting distance to ground and stratification)
[1334]	1265	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]	1266	!
	1267	!-- Limit dissipation rate according to Seifert, Nuijens and
	1268	!-- Stevens (2010)
[1361]	1269	dissipation = MIN( 0.06_wp, diss(k,j,i) )
[1065]	1270	!
	1271	!-- Compute Taylor-microscale Reynolds number:
[1361]	1272	re_lambda = 6.0_wp / 11.0_wp * &
	1273	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	1274	SQRT( 15.0_wp / kin_vis_air ) * &
	1275	dissipation**( 1.0_wp / 6.0_wp )
[1065]	1276	!
	1277	!-- The factor of 1.0E4 is needed to convert the dissipation rate
	1278	!-- from m2 s-3 to cm2 s-3.
[1361]	1279	k_au = k_au * ( 1.0_wp + &
	1280	dissipation * 1.0E4_wp * &
	1281	( re_lambda * 1.0E-3_wp )*0.25_wp &
	1282	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
	1283	sigma_cc )**2 &
	1284	) + beta_cc &
	1285	) &
	1286	)
[1065]	1287	ENDIF
	1288	!
[1012]	1289	!-- Autoconversion rate (Seifert and Beheng, 2006):
[1361]	1290	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	1291	( nu_c + 1.0_wp )*2 qc_1d(k)*2 xc*2 &
	1292	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
[1115]	1293	rho_surface
	1294	autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]	1295
[1115]	1296	qr_1d(k) = qr_1d(k) + autocon * dt_micro
	1297	qc_1d(k) = qc_1d(k) - autocon * dt_micro
	1298	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
	1299
[1005]	1300	ENDIF
[1000]	1301
	1302	ENDDO
	1303
[1005]	1304	END SUBROUTINE autoconversion_ij
	1305
[1106]	1306
[1682]	1307	!------------------------------------------------------------------------------!
	1308	! Description:
	1309	! ------------
	1310	!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
	1311	!------------------------------------------------------------------------------!
[1005]	1312	SUBROUTINE accretion_ij( i, j )
	1313
[1320]	1314	USE arrays_3d, &
	1315	ONLY: diss, qc_1d, qr_1d
[1115]	1316
[1320]	1317	USE cloud_parameters, &
	1318	ONLY: eps_sb, hyrho, k_cr0
	1319
	1320	USE control_parameters, &
	1321	ONLY: dt_micro, rho_surface, turbulence
	1322
	1323	USE indices, &
	1324	ONLY: nzb, nzb_s_inner, nzt
	1325
	1326	USE kinds
	1327
[1005]	1328	IMPLICIT NONE
	1329
[1682]	1330	INTEGER(iwp) :: i !<
	1331	INTEGER(iwp) :: j !<
	1332	INTEGER(iwp) :: k !<
[1005]	1333
[1682]	1334	REAL(wp) :: accr !<
	1335	REAL(wp) :: k_cr !<
	1336	REAL(wp) :: phi_ac !<
	1337	REAL(wp) :: tau_cloud !<
	1338	REAL(wp) :: xc !<
[1320]	1339
[1115]	1340	DO k = nzb_s_inner(j,i)+1, nzt
	1341	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN
[1012]	1342	!
[1048]	1343	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	1344	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
[1012]	1345	!
	1346	!-- Universal function for accretion process
[1048]	1347	!-- (Seifert and Beheng, 2001):
[1361]	1348	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[1012]	1349	!
[1065]	1350	!-- Parameterized turbulence effects on autoconversion (Seifert,
	1351	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]	1352	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1065]	1353	IF ( turbulence ) THEN
[1361]	1354	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	1355	MIN( 600.0_wp, &
	1356	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	1357	)
[1065]	1358	ELSE
	1359	k_cr = k_cr0
	1360	ENDIF
	1361	!
[1012]	1362	!-- Accretion rate (Seifert and Beheng, 2006):
[1361]	1363	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
[1115]	1364	accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]	1365
[1115]	1366	qr_1d(k) = qr_1d(k) + accr * dt_micro
	1367	qc_1d(k) = qc_1d(k) - accr * dt_micro
	1368
[1005]	1369	ENDIF
[1106]	1370
[1005]	1371	ENDDO
	1372
[1000]	1373	END SUBROUTINE accretion_ij
	1374
[1005]	1375
[1682]	1376	!------------------------------------------------------------------------------!
	1377	! Description:
	1378	! ------------
	1379	!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
	1380	!------------------------------------------------------------------------------!
[1005]	1381	SUBROUTINE selfcollection_breakup_ij( i, j )
	1382
[1320]	1383	USE arrays_3d, &
	1384	ONLY: nr_1d, qr_1d
	1385
	1386	USE cloud_parameters, &
	1387	ONLY: dpirho_l, eps_sb, hyrho, k_br, k_rr
	1388
	1389	USE control_parameters, &
	1390	ONLY: dt_micro, rho_surface
	1391
	1392	USE indices, &
	1393	ONLY: nzb, nzb_s_inner, nzt
	1394
	1395	USE kinds
[1005]	1396
	1397	IMPLICIT NONE
	1398
[1682]	1399	INTEGER(iwp) :: i !<
	1400	INTEGER(iwp) :: j !<
	1401	INTEGER(iwp) :: k !<
[1005]	1402
[1682]	1403	REAL(wp) :: breakup !<
	1404	REAL(wp) :: dr !<
	1405	REAL(wp) :: phi_br !<
	1406	REAL(wp) :: selfcoll !<
[1320]	1407
[1115]	1408	DO k = nzb_s_inner(j,i)+1, nzt
	1409	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	1410	!
[1115]	1411	!-- Selfcollection rate (Seifert and Beheng, 2001):
[1361]	1412	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
[1012]	1413	!
[1115]	1414	!-- Weight averaged diameter of rain drops:
[1334]	1415	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	1416	!
[1048]	1417	!-- Collisional breakup rate (Seifert, 2008):
[1353]	1418	IF ( dr >= 0.3E-3_wp ) THEN
	1419	phi_br = k_br * ( dr - 1.1E-3_wp )
	1420	breakup = selfcoll * ( phi_br + 1.0_wp )
[1005]	1421	ELSE
[1353]	1422	breakup = 0.0_wp
[1005]	1423	ENDIF
[1048]	1424
[1115]	1425	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
	1426	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
[1106]	1427
[1005]	1428	ENDIF
	1429	ENDDO
	1430
	1431	END SUBROUTINE selfcollection_breakup_ij
	1432
[1106]	1433
[1682]	1434	!------------------------------------------------------------------------------!
	1435	! Description:
	1436	! ------------
	1437	!> Evaporation of precipitable water. Condensation is neglected for
	1438	!> precipitable water. Call for grid point i,j
	1439	!------------------------------------------------------------------------------!
[1012]	1440	SUBROUTINE evaporation_rain_ij( i, j )
	1441
[1320]	1442	USE arrays_3d, &
	1443	ONLY: hyp, nr_1d, pt_1d, q_1d, qc_1d, qr_1d
[1048]	1444
[1320]	1445	USE cloud_parameters, &
	1446	ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
	1447	dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r, &
	1448	l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l, &
	1449	t_d_pt, ventilation_effect
	1450
	1451	USE constants, &
	1452	ONLY: pi
	1453
	1454	USE control_parameters, &
	1455	ONLY: dt_micro
	1456
	1457	USE indices, &
	1458	ONLY: nzb, nzb_s_inner, nzt
	1459
	1460	USE kinds
	1461
[1012]	1462	IMPLICIT NONE
	1463
[1682]	1464	INTEGER(iwp) :: i !<
	1465	INTEGER(iwp) :: j !<
	1466	INTEGER(iwp) :: k !<
[1012]	1467
[1682]	1468	REAL(wp) :: alpha !<
	1469	REAL(wp) :: dr !<
	1470	REAL(wp) :: e_s !<
	1471	REAL(wp) :: evap !<
	1472	REAL(wp) :: evap_nr !<
	1473	REAL(wp) :: f_vent !<
	1474	REAL(wp) :: g_evap !<
	1475	REAL(wp) :: lambda_r !<
	1476	REAL(wp) :: mu_r !<
	1477	REAL(wp) :: mu_r_2 !<
	1478	REAL(wp) :: mu_r_5d2 !<
	1479	REAL(wp) :: nr_0 !<
	1480	REAL(wp) :: q_s !<
	1481	REAL(wp) :: sat !<
	1482	REAL(wp) :: t_l !<
	1483	REAL(wp) :: temp !<
	1484	REAL(wp) :: xr !<
[1320]	1485
[1115]	1486	DO k = nzb_s_inner(j,i)+1, nzt
	1487	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	1488	!
	1489	!-- Actual liquid water temperature:
[1115]	1490	t_l = t_d_pt(k) * pt_1d(k)
[1012]	1491	!
	1492	!-- Saturation vapor pressure at t_l:
[1361]	1493	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	1494	( t_l - 35.86_wp ) &
	1495	)
[1012]	1496	!
	1497	!-- Computation of saturation humidity:
[1361]	1498	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
[1353]	1499	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
[1361]	1500	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]	1501	!
[1106]	1502	!-- Supersaturation:
[1361]	1503	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
[1012]	1504	!
[1361]	1505	!-- Evaporation needs only to be calculated in subsaturated regions
	1506	IF ( sat < 0.0_wp ) THEN
[1012]	1507	!
[1361]	1508	!-- Actual temperature:
	1509	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
	1510
	1511	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
	1512	( thermal_conductivity_l * temp ) + &
	1513	r_v * temp / ( diff_coeff_l * e_s ) &
	1514	)
[1012]	1515	!
[1361]	1516	!-- Mean weight of rain drops
	1517	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1115]	1518	!
[1361]	1519	!-- Weight averaged diameter of rain drops:
	1520	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	1521	!
[1361]	1522	!-- Compute ventilation factor and intercept parameter
	1523	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	1524	IF ( ventilation_effect ) THEN
[1115]	1525	!
[1361]	1526	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	1527	!-- Stevens and Seifert, 2008):
	1528	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
	1529	!
	1530	!-- Slope parameter of gamma distribution (Seifert, 2008):
	1531	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1532	( mu_r + 1.0_wp ) &
	1533	)**( 1.0_wp / 3.0_wp ) / dr
[1115]	1534
[1361]	1535	mu_r_2 = mu_r + 2.0_wp
	1536	mu_r_5d2 = mu_r + 2.5_wp
	1537
	1538	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
	1539	b_vent * schmidt_p_1d3 * &
	1540	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
	1541	lambda_r*( -mu_r_5d2 ) &
	1542	( 1.0_wp - &
	1543	0.5_wp * ( b_term / a_term ) * &
	1544	( lambda_r / ( c_term + lambda_r ) &
	1545	)**mu_r_5d2 - &
	1546	0.125_wp * ( b_term / a_term )*2 &
	1547	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	1548	)**mu_r_5d2 - &
	1549	0.0625_wp * ( b_term / a_term )*3 &
	1550	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	1551	)**mu_r_5d2 - &
	1552	0.0390625_wp * ( b_term / a_term )*4 &
	1553	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	1554	)**mu_r_5d2 &
	1555	)
	1556
	1557	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
	1558	gamm( mu_r + 1.0_wp )
	1559	ELSE
	1560	f_vent = 1.0_wp
	1561	nr_0 = nr_1d(k) * dr
	1562	ENDIF
[1012]	1563	!
[1361]	1564	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
	1565	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
	1566	evap = MAX( evap, -qr_1d(k) / dt_micro )
	1567	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	1568	-nr_1d(k) / dt_micro )
[1106]	1569
[1361]	1570	qr_1d(k) = qr_1d(k) + evap * dt_micro
	1571	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
[1115]	1572
[1361]	1573	ENDIF
[1012]	1574	ENDIF
[1106]	1575
[1012]	1576	ENDDO
	1577
	1578	END SUBROUTINE evaporation_rain_ij
	1579
[1106]	1580
[1682]	1581	!------------------------------------------------------------------------------!
	1582	! Description:
	1583	! ------------
	1584	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
	1585	!> Call for grid point i,j
	1586	!------------------------------------------------------------------------------!
[1012]	1587	SUBROUTINE sedimentation_cloud_ij( i, j )
	1588
[1320]	1589	USE arrays_3d, &
	1590	ONLY: ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
	1591
	1592	USE cloud_parameters, &
[1361]	1593	ONLY: eps_sb, hyrho, l_d_cp, pt_d_t, sed_qc_const
[1320]	1594
	1595	USE constants, &
	1596	ONLY: pi
	1597
	1598	USE control_parameters, &
	1599	ONLY: dt_do2d_xy, dt_micro, intermediate_timestep_count
	1600
	1601	USE indices, &
	1602	ONLY: nzb, nzb_s_inner, nzt
	1603
	1604	USE kinds
[1012]	1605
	1606	IMPLICIT NONE
	1607
[1682]	1608	INTEGER(iwp) :: i !<
	1609	INTEGER(iwp) :: j !<
	1610	INTEGER(iwp) :: k !<
[1106]	1611
[1682]	1612	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1115]	1613
[1353]	1614	sed_qc(nzt+1) = 0.0_wp
[1012]	1615
[1115]	1616	DO k = nzt, nzb_s_inner(j,i)+1, -1
	1617	IF ( qc_1d(k) > eps_sb ) THEN
[1361]	1618	sed_qc(k) = sed_qc_const * nc_1d(k)*( -2.0_wp / 3.0_wp ) &
	1619	( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
[1115]	1620	ELSE
[1353]	1621	sed_qc(k) = 0.0_wp
[1012]	1622	ENDIF
[1115]	1623
[1361]	1624	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
	1625	dt_micro + sed_qc(k+1) &
	1626	)
[1115]	1627
[1361]	1628	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[1115]	1629	hyrho(k) * dt_micro
[1361]	1630	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[1115]	1631	hyrho(k) * dt_micro
[1361]	1632	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[1115]	1633	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
	1634
[1012]	1635	ENDDO
	1636
	1637	END SUBROUTINE sedimentation_cloud_ij
	1638
[1106]	1639
[1682]	1640	!------------------------------------------------------------------------------!
	1641	! Description:
	1642	! ------------
	1643	!> Computation of sedimentation flux. Implementation according to Stevens
	1644	!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
	1645	!------------------------------------------------------------------------------!
[1012]	1646	SUBROUTINE sedimentation_rain_ij( i, j )
	1647
[1320]	1648	USE arrays_3d, &
	1649	ONLY: ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
	1650
	1651	USE cloud_parameters, &
	1652	ONLY: a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho, &
	1653	limiter_sedimentation, l_d_cp, precipitation_amount, prr, &
	1654	pt_d_t, stp
	1655
	1656	USE control_parameters, &
[1361]	1657	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro, &
	1658	dt_3d, intermediate_timestep_count, &
[1320]	1659	intermediate_timestep_count_max, &
	1660	precipitation_amount_interval, time_do2d_xy
	1661
	1662	USE indices, &
	1663	ONLY: nzb, nzb_s_inner, nzt
	1664
	1665	USE kinds
	1666
	1667	USE statistics, &
	1668	ONLY: weight_substep
[1012]	1669
	1670	IMPLICIT NONE
	1671
[1682]	1672	INTEGER(iwp) :: i !<
	1673	INTEGER(iwp) :: j !<
	1674	INTEGER(iwp) :: k !<
	1675	INTEGER(iwp) :: k_run !<
[1012]	1676
[1682]	1677	REAL(wp) :: c_run !<
	1678	REAL(wp) :: d_max !<
	1679	REAL(wp) :: d_mean !<
	1680	REAL(wp) :: d_min !<
	1681	REAL(wp) :: dr !<
	1682	REAL(wp) :: dt_sedi !<
	1683	REAL(wp) :: flux !<
	1684	REAL(wp) :: lambda_r !<
	1685	REAL(wp) :: mu_r !<
	1686	REAL(wp) :: z_run !<
[1320]	1687
[1682]	1688	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	1689	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	1690	REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !<
	1691	REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !<
	1692	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	1693	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	1694	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	1695	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	1696	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	1697	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1320]	1698
[1353]	1699	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
[1012]	1700	!
[1065]	1701	!-- Compute velocities
	1702	DO k = nzb_s_inner(j,i)+1, nzt
[1115]	1703	IF ( qr_1d(k) > eps_sb ) THEN
	1704	!
	1705	!-- Weight averaged diameter of rain drops:
[1334]	1706	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	1707	!
	1708	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	1709	!-- Stevens and Seifert, 2008):
[1353]	1710	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]	1711	!
	1712	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1361]	1713	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1714	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]	1715
[1361]	1716	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1717	a_term - b_term * ( 1.0_wp + &
	1718	c_term / lambda_r )*( -1.0_wp &
	1719	( mu_r + 1.0_wp ) ) &
	1720	) &
	1721	)
	1722	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1723	a_term - b_term * ( 1.0_wp + &
	1724	c_term / lambda_r )*( -1.0_wp &
	1725	( mu_r + 4.0_wp ) ) &
	1726	) &
	1727	)
[1065]	1728	ELSE
[1353]	1729	w_nr(k) = 0.0_wp
	1730	w_qr(k) = 0.0_wp
[1065]	1731	ENDIF
	1732	ENDDO
[1048]	1733	!
[1065]	1734	!-- Adjust boundary values
[1115]	1735	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
	1736	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
[1353]	1737	w_nr(nzt+1) = 0.0_wp
	1738	w_qr(nzt+1) = 0.0_wp
[1065]	1739	!
	1740	!-- Compute Courant number
[1115]	1741	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	1742	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[1115]	1743	dt_micro * ddzu(k)
[1361]	1744	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[1115]	1745	dt_micro * ddzu(k)
	1746	ENDDO
[1065]	1747	!
	1748	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	1749	IF ( limiter_sedimentation ) THEN
	1750
[1115]	1751	DO k = nzb_s_inner(j,i)+1, nzt
[1646]	1752	d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
[1115]	1753	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
	1754	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]	1755
[1361]	1756	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1757	2.0_wp * d_max, &
	1758	ABS( d_mean ) )
[1065]	1759
[1646]	1760	d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
[1115]	1761	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
	1762	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]	1763
[1361]	1764	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1765	2.0_wp * d_max, &
	1766	ABS( d_mean ) )
[1022]	1767	ENDDO
[1048]	1768
[1065]	1769	ELSE
[1106]	1770
[1353]	1771	nr_slope = 0.0_wp
	1772	qr_slope = 0.0_wp
[1106]	1773
[1065]	1774	ENDIF
[1115]	1775
[1353]	1776	sed_nr(nzt+1) = 0.0_wp
	1777	sed_qr(nzt+1) = 0.0_wp
[1065]	1778	!
	1779	!-- Compute sedimentation flux
[1115]	1780	DO k = nzt, nzb_s_inner(j,i)+1, -1
[1065]	1781	!
	1782	!-- Sum up all rain drop number densities which contribute to the flux
	1783	!-- through k-1/2
[1353]	1784	flux = 0.0_wp
	1785	z_run = 0.0_wp ! height above z(k)
[1065]	1786	k_run = k
[1346]	1787	c_run = MIN( 1.0_wp, c_nr(k) )
[1353]	1788	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1361]	1789	flux = flux + hyrho(k_run) * &
	1790	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
[1353]	1791	0.5_wp ) * c_run * dzu(k_run)
[1065]	1792	z_run = z_run + dzu(k_run)
	1793	k_run = k_run + 1
[1346]	1794	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
[1022]	1795	ENDDO
	1796	!
[1065]	1797	!-- It is not allowed to sediment more rain drop number density than
	1798	!-- available
[1361]	1799	flux = MIN( flux, &
[1115]	1800	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]	1801
[1115]	1802	sed_nr(k) = flux / dt_micro
[1361]	1803	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
	1804	hyrho(k) * dt_micro
[1065]	1805	!
	1806	!-- Sum up all rain water content which contributes to the flux
	1807	!-- through k-1/2
[1353]	1808	flux = 0.0_wp
	1809	z_run = 0.0_wp ! height above z(k)
[1065]	1810	k_run = k
[1346]	1811	c_run = MIN( 1.0_wp, c_qr(k) )
[1106]	1812
[1361]	1813	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1106]	1814
[1361]	1815	flux = flux + hyrho(k_run) * &
	1816	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
[1353]	1817	0.5_wp ) * c_run * dzu(k_run)
[1065]	1818	z_run = z_run + dzu(k_run)
	1819	k_run = k_run + 1
[1346]	1820	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
[1106]	1821
[1065]	1822	ENDDO
	1823	!
	1824	!-- It is not allowed to sediment more rain water content than available
[1361]	1825	flux = MIN( flux, &
[1115]	1826	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]	1827
[1115]	1828	sed_qr(k) = flux / dt_micro
	1829
[1361]	1830	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[1115]	1831	hyrho(k) * dt_micro
[1361]	1832	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[1115]	1833	hyrho(k) * dt_micro
[1361]	1834	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[1115]	1835	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
[1065]	1836	!
	1837	!-- Compute the rain rate
[1361]	1838	IF ( call_microphysics_at_all_substeps ) THEN
	1839	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
	1840	weight_substep(intermediate_timestep_count)
	1841	ELSE
	1842	prr(k,j,i) = sed_qr(k) / hyrho(k)
	1843	ENDIF
	1844
[1065]	1845	ENDDO
[1115]	1846
[1065]	1847	!
[1048]	1848	!-- Precipitation amount
[1361]	1849	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
	1850	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
	1851	precipitation_amount_interval ) THEN
[1012]	1852
[1361]	1853	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	1854	prr(nzb_s_inner(j,i)+1,j,i) * &
[1115]	1855	hyrho(nzb_s_inner(j,i)+1) * dt_3d
[1048]	1856	ENDIF
	1857
[1012]	1858	END SUBROUTINE sedimentation_rain_ij
	1859
[1682]	1860
[1361]	1861	!------------------------------------------------------------------------------!
[1682]	1862	! Description:
	1863	! ------------
	1864	!> This function computes the gamma function (Press et al., 1992).
	1865	!> The gamma function is needed for the calculation of the evaporation
	1866	!> of rain drops.
[1361]	1867	!------------------------------------------------------------------------------!
[1012]	1868	FUNCTION gamm( xx )
[1048]	1869
[1320]	1870	USE cloud_parameters, &
	1871	ONLY: cof, stp
	1872
	1873	USE kinds
	1874
[1012]	1875	IMPLICIT NONE
[1106]	1876
[1682]	1877	INTEGER(iwp) :: j !<
[1320]	1878
[1682]	1879	REAL(wp) :: gamm !<
	1880	REAL(wp) :: ser !<
	1881	REAL(wp) :: tmp !<
	1882	REAL(wp) :: x_gamm !<
	1883	REAL(wp) :: xx !<
	1884	REAL(wp) :: y_gamm !<
[1320]	1885
[1012]	1886	x_gamm = xx
	1887	y_gamm = x_gamm
[1353]	1888	tmp = x_gamm + 5.5_wp
	1889	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
[1334]	1890	ser = 1.000000000190015_wp
[1106]	1891
	1892	DO j = 1, 6
[1353]	1893	y_gamm = y_gamm + 1.0_wp
[1012]	1894	ser = ser + cof( j ) / y_gamm
[1106]	1895	ENDDO
	1896
[1012]	1897	!
	1898	!-- Until this point the algorithm computes the logarithm of the gamma
	1899	!-- function. Hence, the exponential function is used.
	1900	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
	1901	gamm = EXP( tmp ) * stp * ser / x_gamm
[1106]	1902
[1012]	1903	RETURN
	1904
	1905	END FUNCTION gamm
	1906
	1907	END MODULE microphysics_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |