Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |