Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_diffusion_em.F @ 3094

Last change on this file since 3094 was 2761, checked in by aslmd, 2 years ago
Applied planetary adaptation changes to WRFV3. job done previously by LMD_LES_MARS_install. Moved Registry.EM.newphys to put it simply in Registry.EM
File size: 193.3 KB

Rev	Line
[17]	1	! WRF:MODEL_LAYER:PHYSICS
	2
	3	MODULE module_diffusion_em
	4
	5	USE module_configure
	6	USE module_bc
	7	USE module_state_description
	8	USE module_big_step_utilities_em
	9	USE module_model_constants
	10	USE module_wrf_error
	11
	12	CONTAINS
	13
	14	!=======================================================================
	15	!=======================================================================
	16
	17	SUBROUTINE cal_deform_and_div( config_flags, u, v, w, div, &
	18	defor11, defor22, defor33, &
	19	defor12, defor13, defor23, &
	20	u_base, v_base, msfux, msfuy, &
	21	msfvx, msfvy, msftx, msfty, &
	22	rdx, rdy, dn, dnw, rdz, rdzw, &
	23	fnm, fnp, cf1, cf2, cf3, zx, zy, &
	24	ids, ide, jds, jde, kds, kde, &
	25	ims, ime, jms, jme, kms, kme, &
	26	its, ite, jts, jte, kts, kte )
	27
	28	! History: Sep 2003 Changes by Jason Knievel and George Bryan, NCAR
	29	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	30	! ... ...
	31
	32	! Purpose: This routine calculates deformation and 3-d divergence.
	33
	34	! References: Klemp and Wilhelmson (JAS 1978)
	35	! Chen and Dudhia (NCAR WRF physics report 2000)
	36
	37	!-----------------------------------------------------------------------
	38	! Comments 10-MAR-05
	39	! Equations 13a-f, Chen and Dudhia 2000, Appendix A:
	40	! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
	41	! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
	42	! Eqn 13c: D33=defor33= 2 * partial dw/dz [SIMPLER FORM]
	43	! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
	44	! partial dpsi/dx * partial dv^/dpsi +
	45	! partial dpsi/dy * partial du^/dpsi)
	46	! Eqn 13e: D13=defor13= m^2 * (partial dw^/dX + partial dpsi/dx * partial dw^/dpsi)
	47	! + partial du/dz [SIMPLER FORM]
	48	! Eqn 13f: D23=defor23= m^2 * (partial dw^/dY + partial dpsi/dy * partial dw^/dpsi)
	49	! + partial dv/dz [SIMPLER FORM]
	50	!-----------------------------------------------------------------------
	51	! Begin declarations.
	52
	53	IMPLICIT NONE
	54
	55	TYPE( grid_config_rec_type ), INTENT( IN ) &
	56	:: config_flags
	57
	58	INTEGER, INTENT( IN ) &
	59	:: ids, ide, jds, jde, kds, kde, &
	60	ims, ime, jms, jme, kms, kme, &
	61	its, ite, jts, jte, kts, kte
	62
	63	REAL, INTENT( IN ) &
	64	:: rdx, rdy, cf1, cf2, cf3
	65
	66	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
	67	:: fnm, fnp, dn, dnw, u_base, v_base
	68
	69	REAL, DIMENSION( ims:ime , jms:jme ), INTENT( IN ) &
	70	:: msfux, msfuy, msfvx, msfvy, msftx, msfty
	71
	72	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	73	:: u, v, w, zx, zy, rdz, rdzw
	74
	75	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	76	:: defor11, defor22, defor33, defor12, defor13, defor23, div
	77
	78	! Local variables.
	79
	80	INTEGER &
	81	:: i, j, k, ktf, ktes1, ktes2, i_start, i_end, j_start, j_end
	82
	83	REAL &
	84	:: tmp, tmpzx, tmpzy, tmpzeta_z, cft1, cft2
	85
	86	REAL, DIMENSION( its:ite, jts:jte ) &
	87	:: mm, zzavg, zeta_zd12
	88
	89	REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) &
	90	:: tmp1, hat, hatavg
	91
	92	! End declarations.
	93	!-----------------------------------------------------------------------
	94
	95	! Comments 10-MAR-2005
	96	! Treat all differentials as 'div-style' [or 'curl-style'],
	97	! i.e., du/dX becomes (in map coordinate space) mxmy d(u/my)/dx,
	98	! NB - all equations referred to here are from Chen and Dudhia 2002, from the
	99	! WRF physics documents web pages:
	100	! http://www.mmm.ucar.edu/wrf/users/docs/wrf-doc-physics.pdf
	101
	102	!=======================================================================
	103	! In the following section, calculate 3-d divergence and the first three
	104	! (defor11, defor22, defor33) of six deformation terms.
	105
	106	ktes1 = kte-1
	107	ktes2 = kte-2
	108
	109	cft2 = - 0.5 * dnw(ktes1) / dn(ktes1)
	110	cft1 = 1.0 - cft2
	111
	112	ktf = MIN( kte, kde-1 )
	113
	114	i_start = its
	115	i_end = MIN( ite, ide-1 )
	116	j_start = jts
	117	j_end = MIN( jte, jde-1 )
	118
	119	! Square the map scale factor.
	120
	121	DO j = j_start, j_end
	122	DO i = i_start, i_end
	123	mm(i,j) = msftx(i,j) * msfty(i,j)
	124	END DO
	125	END DO
	126
	127	!-----------------------------------------------------------------------
	128	! Calculate du/dx.
	129
	130	! Apply a coordinate transformation to zonal velocity, u.
	131
	132	DO j = j_start, j_end
	133	DO k = kts, ktf
	134	DO i = i_start, i_end+1
	135	hat(i,k,j) = u(i,k,j) / msfuy(i,j)
	136	END DO
	137	END DO
	138	END DO
	139
	140	! Average in x and z.
	141
	142	DO j=j_start,j_end
	143	DO k=kts+1,ktf
	144	DO i=i_start,i_end
	145	hatavg(i,k,j) = 0.5 * &
	146	( fnm(k) * ( hat(i,k ,j) + hat(i+1, k,j) ) + &
	147	fnp(k) * ( hat(i,k-1,j) + hat(i+1,k-1,j) ) )
	148	END DO
	149	END DO
	150	END DO
	151
	152	! Extrapolate to top and bottom of domain (to w levels).
	153
	154	DO j = j_start, j_end
	155	DO i = i_start, i_end
	156	hatavg(i,1,j) = 0.5 * ( &
	157	cf1 * hat(i ,1,j) + &
	158	cf2 * hat(i ,2,j) + &
	159	cf3 * hat(i ,3,j) + &
	160	cf1 * hat(i+1,1,j) + &
	161	cf2 * hat(i+1,2,j) + &
	162	cf3 * hat(i+1,3,j) )
	163	hatavg(i,kte,j) = 0.5 * ( &
	164	cft1 * ( hat(i,ktes1,j) + hat(i+1,ktes1,j) ) + &
	165	cft2 * ( hat(i,ktes2,j) + hat(i+1,ktes2,j) ) )
	166	END DO
	167	END DO
	168
	169	! Comments 10-MAR-05
	170	! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
	171	! Below, D11 is set = 2*tmp1
	172	! => tmp1 = m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
	173	! tmpzx = averaged value of dpsi/dx (=zx)
	174
	175	DO j = j_start, j_end
	176	DO k = kts, ktf
	177	DO i = i_start, i_end
	178	tmpzx = 0.25 * ( &
	179	zx(i,k ,j) + zx(i+1,k ,j) + &
	180	zx(i,k+1,j) + zx(i+1,k+1,j) )
	181	tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) tmpzx rdzw(i,k,j)
	182	! tmp1 to here = partial dpsi/dx * partial du^/dpsi:
	183	END DO
	184	END DO
	185	END DO
	186
	187	DO j = j_start, j_end
	188	DO k = kts, ktf
	189	DO i = i_start, i_end
	190	tmp1(i,k,j) = mm(i,j) * ( rdx * ( hat(i+1,k,j) - hat(i,k,j) ) - &
	191	tmp1(i,k,j))
	192	END DO
	193	END DO
	194	END DO
	195
	196	! End calculation of du/dx.
	197	!-----------------------------------------------------------------------
	198
	199	!-----------------------------------------------------------------------
	200	! Calculate defor11 (2*du/dx).
	201	! Comments 10-MAR-05
	202	! Eqn 13a: D11=defor11= 2 m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
	203	! = 2*tmp1
	204
	205	DO j = j_start, j_end
	206	DO k = kts, ktf
	207	DO i = i_start, i_end
	208	defor11(i,k,j) = 2.0 * tmp1(i,k,j)
	209	END DO
	210	END DO
	211	END DO
	212
	213	! End calculation of defor11.
	214	!-----------------------------------------------------------------------
	215
	216	!-----------------------------------------------------------------------
	217	! Calculate zonal divergence (du/dx) and add it to the divergence array.
	218
	219	DO j = j_start, j_end
	220	DO k = kts, ktf
	221	DO i = i_start, i_end
	222	div(i,k,j) = tmp1(i,k,j)
	223	END DO
	224	END DO
	225	END DO
	226
	227	! End calculation of zonal divergence.
	228	!-----------------------------------------------------------------------
	229
	230	!-----------------------------------------------------------------------
	231	! Calculate dv/dy.
	232
	233	! Apply a coordinate transformation to meridional velocity, v.
	234
	235	DO j = j_start, j_end+1
	236	DO k = kts, ktf
	237	DO i = i_start, i_end
	238	! Because msfvx at the poles will be undefined (1./0.), we will have
	239	! trouble. But we are OK since v at the poles is 0., and that takes
	240	! precedence in this case.
	241	IF ((config_flags%polar) .AND. ((j == jds) .OR. (j == jde))) THEN
	242	hat(i,k,j) = 0.
	243	ELSE ! normal code
	244	hat(i,k,j) = v(i,k,j) / msfvx(i,j)
	245	ENDIF
	246	END DO
	247	END DO
	248	END DO
	249
	250	! Account for the slope in y of eta surfaces.
	251
	252	DO j=j_start,j_end
	253	DO k=kts+1,ktf
	254	DO i=i_start,i_end
	255	hatavg(i,k,j) = 0.5 * ( &
	256	fnm(k) * ( hat(i,k ,j) + hat(i,k ,j+1) ) + &
	257	fnp(k) * ( hat(i,k-1,j) + hat(i,k-1,j+1) ) )
	258	END DO
	259	END DO
	260	END DO
	261
	262	! Extrapolate to top and bottom of domain (to w levels).
	263
	264	DO j = j_start, j_end
	265	DO i = i_start, i_end
	266	hatavg(i,1,j) = 0.5 * ( &
	267	cf1 * hat(i,1,j ) + &
	268	cf2 * hat(i,2,j ) + &
	269	cf3 * hat(i,3,j ) + &
	270	cf1 * hat(i,1,j+1) + &
	271	cf2 * hat(i,2,j+1) + &
	272	cf3 * hat(i,3,j+1) )
	273	hatavg(i,kte,j) = 0.5 * ( &
	274	cft1 * ( hat(i,ktes1,j) + hat(i,ktes1,j+1) ) + &
	275	cft2 * ( hat(i,ktes2,j) + hat(i,ktes2,j+1) ) )
	276	END DO
	277	END DO
	278
	279	! Comments 10-MAR-05
	280	! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
	281	! Below, D22 is set = 2*tmp1
	282	! => tmp1 = m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
	283	! tmpzy = averaged value of dpsi/dy (=zy)
	284
	285	DO j = j_start, j_end
	286	DO k = kts, ktf
	287	DO i = i_start, i_end
	288	tmpzy = 0.25 * ( &
	289	zy(i,k ,j) + zy(i,k ,j+1) + &
	290	zy(i,k+1,j) + zy(i,k+1,j+1) )
	291	tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * tmpzy * rdzw(i,k,j)
	292	! tmp1 to here = partial dpsi/dy * partial dv^/dpsi:
	293	END DO
	294	END DO
	295	END DO
	296
	297	DO j = j_start, j_end
	298	DO k = kts, ktf
	299	DO i = i_start, i_end
	300	tmp1(i,k,j) = mm(i,j) * ( &
	301	rdy * ( hat(i,k,j+1) - hat(i,k,j) ) - tmp1(i,k,j) )
	302	END DO
	303	END DO
	304	END DO
	305
	306	! End calculation of dv/dy.
	307	!-----------------------------------------------------------------------
	308
	309	!-----------------------------------------------------------------------
	310	! Calculate defor22 (2*dv/dy).
	311	! Comments 10-MAR-05
	312	! Eqn 13b: D22=defor22= 2 m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
	313	! = 2*tmp1
	314
	315	DO j = j_start, j_end
	316	DO k = kts, ktf
	317	DO i = i_start, i_end
	318	defor22(i,k,j) = 2.0 * tmp1(i,k,j)
	319	END DO
	320	END DO
	321	END DO
	322
	323	! End calculation of defor22.
	324	!-----------------------------------------------------------------------
	325
	326	!-----------------------------------------------------------------------
	327	! Calculate meridional divergence (dv/dy) and add it to the divergence
	328	! array.
	329
	330	DO j = j_start, j_end
	331	DO k = kts, ktf
	332	DO i = i_start, i_end
	333	div(i,k,j) = div(i,k,j) + tmp1(i,k,j)
	334	END DO
	335	END DO
	336	END DO
	337
	338	! End calculation of meridional divergence.
	339	!-----------------------------------------------------------------------
	340
	341	!-----------------------------------------------------------------------
	342	! Comments 10-MAR-05
	343	! Eqn 13c: D33=defor33= 2 * partial dw/dz
	344	! Below, D33 is set = 2*tmp1
	345	! => tmp1 = partial dw/dz
	346
	347	! Calculate dw/dz.
	348
	349	DO j = j_start, j_end
	350	DO k = kts, ktf
	351	DO i = i_start, i_end
	352	tmp1(i,k,j) = ( w(i,k+1,j) - w(i,k,j) ) * rdzw(i,k,j)
	353	END DO
	354	END DO
	355	END DO
	356
	357	! End calculation of dw/dz.
	358	!-----------------------------------------------------------------------
	359
	360	!-----------------------------------------------------------------------
	361	! Calculate defor33 (2*dw/dz).
	362
	363	DO j = j_start, j_end
	364	DO k = kts, ktf
	365	DO i = i_start, i_end
	366	defor33(i,k,j) = 2.0 * tmp1(i,k,j)
	367	END DO
	368	END DO
	369	END DO
	370
	371	! End calculation of defor33.
	372	!-----------------------------------------------------------------------
	373
	374	!-----------------------------------------------------------------------
	375	! Calculate vertical divergence (dw/dz) and add it to the divergence
	376	! array.
	377
	378	DO j = j_start, j_end
	379	DO k = kts, ktf
	380	DO i = i_start, i_end
	381	div(i,k,j) = div(i,k,j) + tmp1(i,k,j)
	382	END DO
	383	END DO
	384	END DO
	385
	386	! End calculation of vertical divergence.
	387	!-----------------------------------------------------------------------
	388
	389	! Three-dimensional divergence is now finished and values are in array
	390	! "div." Also, the first three (defor11, defor22, defor33) of six
	391	! deformation terms are now calculated at pressure points.
	392	!=======================================================================
	393
	394	! Comments 10-MAR-2005
	395	! Treat all differentials as 'div-style' [or 'curl-style'],
	396	! i.e., du/dY becomes (in map coordinate space) mxmy d(u/mx)/dy,
	397	! dv/dX becomes (in map coordinate space) mxmy d(v/my)/dx,
	398	! (see e.g. Haltiner and Williams p. 441)
	399
	400	!=======================================================================
	401	! Calculate the final three deformations (defor12, defor13, defor23) at
	402	! vorticity points.
	403
	404	i_start = its
	405	i_end = ite
	406	j_start = jts
	407	j_end = jte
	408
	409	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	410	config_flags%nested) i_start = MAX( ids+1, its )
	411	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	412	config_flags%nested) i_end = MIN( ide-1, ite )
	413	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	414	config_flags%nested) j_start = MAX( jds+1, jts )
	415	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	416	config_flags%nested) j_end = MIN( jde-1, jte )
	417	IF ( config_flags%periodic_x ) i_start = its
	418	IF ( config_flags%periodic_x ) i_end = ite
	419
	420
	421	!-----------------------------------------------------------------------
	422	! Calculate du/dy.
	423
	424	! First, calculate an average mapscale factor.
	425
	426	! Comments 10-MAR-05
	427	! du/dy => need u map scale factor in x (which is defined at u points)
	428	! averaged over j and j-1
	429	! dv/dx => need v map scale factor in y (which is defined at v points)
	430	! averaged over i and i-1
	431
	432	DO j = j_start, j_end
	433	DO i = i_start, i_end
	434	mm(i,j) = 0.25 * ( msfux(i,j-1) + msfux(i,j) ) * ( msfvy(i-1,j) + msfvy(i,j) )
	435	END DO
	436	END DO
	437
	438	! Apply a coordinate transformation to zonal velocity, u.
	439
	440	DO j =j_start-1, j_end
	441	DO k =kts, ktf
	442	DO i =i_start, i_end
	443	! Fixes to set_physical_bc2/3d for polar boundary conditions
	444	! remove issues with loop over j
	445	hat(i,k,j) = u(i,k,j) / msfux(i,j)
	446	END DO
	447	END DO
	448	END DO
	449
	450	! Average in y and z.
	451
	452	DO j=j_start,j_end
	453	DO k=kts+1,ktf
	454	DO i=i_start,i_end
	455	hatavg(i,k,j) = 0.5 * ( &
	456	fnm(k) * ( hat(i,k ,j-1) + hat(i,k ,j) ) + &
	457	fnp(k) * ( hat(i,k-1,j-1) + hat(i,k-1,j) ) )
	458	END DO
	459	END DO
	460	END DO
	461
	462	! Extrapolate to top and bottom of domain (to w levels).
	463
	464	DO j = j_start, j_end
	465	DO i = i_start, i_end
	466	hatavg(i,1,j) = 0.5 * ( &
	467	cf1 * hat(i,1,j-1) + &
	468	cf2 * hat(i,2,j-1) + &
	469	cf3 * hat(i,3,j-1) + &
	470	cf1 * hat(i,1,j ) + &
	471	cf2 * hat(i,2,j ) + &
	472	cf3 * hat(i,3,j ) )
	473	hatavg(i,kte,j) = 0.5 * ( &
	474	cft1 * ( hat(i,ktes1,j-1) + hat(i,ktes1,j) ) + &
	475	cft2 * ( hat(i,ktes2,j-1) + hat(i,ktes2,j) ) )
	476	END DO
	477	END DO
	478
	479	! tmpzy = averaged value of dpsi/dy (=zy) on vorticity grid
	480	! tmp1 = partial dpsi/dy * partial du^/dpsi
	481	DO j = j_start, j_end
	482	DO k = kts, ktf
	483	DO i = i_start, i_end
	484	tmpzy = 0.25 * ( &
	485	zy(i-1,k ,j) + zy(i,k ,j) + &
	486	zy(i-1,k+1,j) + zy(i,k+1,j) )
	487	tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * &
	488	0.25 * tmpzy * ( rdzw(i,k,j) + rdzw(i-1,k,j) + &
	489	rdzw(i-1,k,j-1) + rdzw(i,k,j-1) )
	490	END DO
	491	END DO
	492	END DO
	493
	494	! End calculation of du/dy.
	495	!----------------------------------------------------------------------
	496
	497	!-----------------------------------------------------------------------
	498	! Add the first term to defor12 (du/dy+dv/dx) at vorticity points.
	499
	500	! Comments 10-MAR-05
	501	! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
	502	! partial dpsi/dx * partial dv^/dpsi +
	503	! partial dpsi/dy * partial du^/dpsi)
	504	! Here deal with m^2 * (partial du^/dY + partial dpsi/dy * partial du^/dpsi)
	505	! Still need to add v^ terms:
	506	! m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi)
	507
	508	DO j = j_start, j_end
	509	DO k = kts, ktf
	510	DO i = i_start, i_end
	511	defor12(i,k,j) = mm(i,j) * ( &
	512	rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) )
	513	END DO
	514	END DO
	515	END DO
	516
	517	! End addition of the first term to defor12.
	518	!-----------------------------------------------------------------------
	519
	520	!-----------------------------------------------------------------------
	521	! Calculate dv/dx.
	522
	523	! Apply a coordinate transformation to meridional velocity, v.
	524
	525	DO j = j_start, j_end
	526	DO k = kts, ktf
	527	DO i = i_start-1, i_end
	528	hat(i,k,j) = v(i,k,j) / msfvy(i,j)
	529	END DO
	530	END DO
	531	END DO
	532
	533	! Account for the slope in x of eta surfaces.
	534
	535	DO j = j_start, j_end
	536	DO k = kts+1, ktf
	537	DO i = i_start, i_end
	538	hatavg(i,k,j) = 0.5 * ( &
	539	fnm(k) * ( hat(i-1,k ,j) + hat(i,k ,j) ) + &
	540	fnp(k) * ( hat(i-1,k-1,j) + hat(i,k-1,j) ) )
	541	END DO
	542	END DO
	543	END DO
	544
	545	! Extrapolate to top and bottom of domain (to w levels).
	546
	547	DO j = j_start, j_end
	548	DO i = i_start, i_end
	549	hatavg(i,1,j) = 0.5 * ( &
	550	cf1 * hat(i-1,1,j) + &
	551	cf2 * hat(i-1,2,j) + &
	552	cf3 * hat(i-1,3,j) + &
	553	cf1 * hat(i ,1,j) + &
	554	cf2 * hat(i ,2,j) + &
	555	cf3 * hat(i ,3,j) )
	556	hatavg(i,kte,j) = 0.5 * ( &
	557	cft1 * ( hat(i,ktes1,j) + hat(i-1,ktes1,j) ) + &
	558	cft2 * ( hat(i,ktes2,j) + hat(i-1,ktes2,j) ) )
	559	END DO
	560	END DO
	561
	562	! Fixes to set_physical_bc2/3d have made any check for polar B.C.'s
	563	! unnecessary in this place. zx, rdzw, and hatavg are all defined
	564	! in places they need to be and the values at the poles are replications
	565	! of the values one grid point in, so the averaging over j and j-1 works
	566	! to act as just using the value at j or j-1 (with out extra code).
	567	!
	568	! tmpzx = averaged value of dpsi/dx (=zx) on vorticity grid
	569	! tmp1 = partial dpsi/dx * partial dv^/dpsi
	570	DO j = j_start, j_end
	571	DO k = kts, ktf
	572	DO i = i_start, i_end
	573	tmpzx = 0.25 * ( &
	574	zx(i,k ,j-1) + zx(i,k ,j) + &
	575	zx(i,k+1,j-1) + zx(i,k+1,j) )
	576	tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * &
	577	0.25 * tmpzx * ( rdzw(i,k,j) + rdzw(i,k,j-1) + &
	578	rdzw(i-1,k,j-1) + rdzw(i-1,k,j) )
	579	END DO
	580	END DO
	581	END DO
	582
	583	! End calculation of dv/dx.
	584	!-----------------------------------------------------------------------
	585
	586	!-----------------------------------------------------------------------
	587	! Add the second term to defor12 (du/dy+dv/dx) at vorticity points.
	588
	589	! Comments 10-MAR-05
	590	! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
	591	! partial dpsi/dx * partial dv^/dpsi +
	592	! partial dpsi/dy * partial du^/dpsi)
	593	! Here adding v^ terms:
	594	! m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi)
	595
	596	DO j = j_start, j_end
	597	DO k = kts, ktf
	598	DO i = i_start, i_end
	599	defor12(i,k,j) = defor12(i,k,j) + &
	600	mm(i,j) * ( &
	601	rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
	602	END DO
	603	END DO
	604	END DO
	605
	606	! End addition of the second term to defor12.
	607	!-----------------------------------------------------------------------
	608
	609	!-----------------------------------------------------------------------
	610	! Update the boundary for defor12 (might need to change later).
	611
	612	IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1 ) THEN
	613	DO j = jts, jte
	614	DO k = kts, kte
	615	defor12(ids,k,j) = defor12(ids+1,k,j)
	616	END DO
	617	END DO
	618	END IF
	619
	620	IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
	621	DO k = kts, kte
	622	DO i = its, ite
	623	defor12(i,k,jds) = defor12(i,k,jds+1)
	624	END DO
	625	END DO
	626	END IF
	627
	628	IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
	629	DO j = jts, jte
	630	DO k = kts, kte
	631	defor12(ide,k,j) = defor12(ide-1,k,j)
	632	END DO
	633	END DO
	634	END IF
	635
	636	IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
	637	DO k = kts, kte
	638	DO i = its, ite
	639	defor12(i,k,jde) = defor12(i,k,jde-1)
	640	END DO
	641	END DO
	642	END IF
	643
	644	! End update of boundary for defor12.
	645	!-----------------------------------------------------------------------
	646
	647	! Comments 10-MAR-05
	648	! Further deformation terms not needed for 2-dimensional Smagorinsky diffusion,
	649	! so those terms have not been dealt with yet.
	650	! A "y" has simply been added to all map scale factors to allow the model to
	651	! compile without errors.
	652
	653	!-----------------------------------------------------------------------
	654	! Calculate dw/dx.
	655
	656	i_start = its
	657	i_end = MIN( ite, ide-1 )
	658	j_start = jts
	659	j_end = MIN( jte, jde-1 )
	660
	661	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	662	config_flags%nested) i_start = MAX( ids+1, its )
	663	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	664	config_flags%nested) j_start = MAX( jds+1, jts )
	665
	666	IF ( config_flags%periodic_x ) i_start = its
	667	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide )
	668	IF ( config_flags%periodic_y ) j_end = MIN( jte, jde )
	669
	670	! Square the mapscale factor.
	671
	672	DO j = jts, jte
	673	DO i = its, ite
	674	mm(i,j) = msfux(i,j) * msfuy(i,j)
	675	END DO
	676	END DO
	677
	678	! Apply a coordinate transformation to vertical velocity, w. This is for both
	679	! defor13 and defor23.
	680
	681	DO j = j_start, j_end
	682	DO k = kts, kte
	683	DO i = i_start, i_end
	684	hat(i,k,j) = w(i,k,j) / msfty(i,j)
	685	END DO
	686	END DO
	687	END DO
	688
	689	i = i_start-1
	690	DO j = j_start, MIN( jte, jde-1 )
	691	DO k = kts, kte
	692	hat(i,k,j) = w(i,k,j) / msfty(i,j)
	693	END DO
	694	END DO
	695
	696	j = j_start-1
	697	DO k = kts, kte
	698	DO i = i_start, MIN( ite, ide-1 )
	699	hat(i,k,j) = w(i,k,j) / msfty(i,j)
	700	END DO
	701	END DO
	702
	703	! QUESTION: What is this for?
	704
	705	DO j = j_start, j_end
	706	DO k = kts, ktf
	707	DO i = i_start, i_end
	708	hatavg(i,k,j) = 0.25 * ( &
	709	hat(i ,k ,j) + &
	710	hat(i ,k+1,j) + &
	711	hat(i-1,k ,j) + &
	712	hat(i-1,k+1,j) )
	713	END DO
	714	END DO
	715	END DO
	716
	717	! Calculate dw/dx.
	718
	719	DO j = j_start, j_end
	720	DO k = kts+1, ktf
	721	DO i = i_start, i_end
	722	tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zx(i,k,j) * &
	723	0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
	724	END DO
	725	END DO
	726	END DO
	727
	728	! End calculation of dw/dx.
	729	!-----------------------------------------------------------------------
	730
	731	!-----------------------------------------------------------------------
	732	! Add the first term (dw/dx) to defor13 (dw/dx+du/dz) at vorticity
	733	! points.
	734
	735	DO j = j_start, j_end
	736	DO k = kts+1, ktf
	737	DO i = i_start, i_end
	738	defor13(i,k,j) = mm(i,j) * ( &
	739	rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
	740	END DO
	741	END DO
	742	END DO
	743
	744	DO j = j_start, j_end
	745	DO i = i_start, i_end
	746	defor13(i,kts,j ) = 0.0
	747	defor13(i,ktf+1,j) = 0.0
	748	END DO
	749	END DO
	750
	751	! End addition of the first term to defor13.
	752	!-----------------------------------------------------------------------
	753
	754	!-----------------------------------------------------------------------
	755	! Calculate du/dz.
	756
	757	IF ( config_flags%mix_full_fields ) THEN
	758
	759	DO j = j_start, j_end
	760	DO k = kts+1, ktf
	761	DO i = i_start, i_end
	762	tmp1(i,k,j) = ( u(i,k,j) - u(i,k-1,j) ) * &
	763	0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
	764	END DO
	765	END DO
	766	END DO
	767
	768	ELSE
	769
	770	DO j = j_start, j_end
	771	DO k = kts+1, ktf
	772	DO i = i_start, i_end
	773	tmp1(i,k,j) = ( u(i,k,j) - u_base(k) - u(i,k-1,j) + u_base(k-1) ) * &
	774	0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
	775	END DO
	776	END DO
	777	END DO
	778
	779	END IF
	780
	781	!-----------------------------------------------------------------------
	782	! Add the second term (du/dz) to defor13 (dw/dx+du/dz) at vorticity
	783	! points.
	784
	785	DO j = j_start, j_end
	786	DO k = kts+1, ktf
	787	DO i = i_start, i_end
	788	defor13(i,k,j) = defor13(i,k,j) + tmp1(i,k,j)
	789	END DO
	790	END DO
	791	END DO
	792
	793	! End addition of the second term to defor13.
	794	!-----------------------------------------------------------------------
	795
	796	!-----------------------------------------------------------------------
	797	! Calculate dw/dy.
	798
	799	i_start = its
	800	i_end = MIN( ite, ide-1 )
	801	j_start = jts
	802	j_end = MIN( jte, jde-1 )
	803
	804	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	805	config_flags%nested) i_start = MAX( ids+1, its )
	806	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	807	config_flags%nested) j_start = MAX( jds+1, jts )
	808	IF ( config_flags%periodic_y ) j_end = MIN( jte, jde )
	809	IF ( config_flags%periodic_x ) i_start = its
	810	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	811
	812	! Square mapscale factor.
	813
	814	DO j = jts, jte
	815	DO i = its, ite
	816	mm(i,j) = msfvx(i,j) * msfvy(i,j)
	817	END DO
	818	END DO
	819
	820	! Apply a coordinate transformation to vertical velocity, w. Added by CW 7/19/07
	821
	822	DO j = j_start, j_end
	823	DO k = kts, kte
	824	DO i = i_start, i_end
	825	hat(i,k,j) = w(i,k,j) / msftx(i,j)
	826	END DO
	827	END DO
	828	END DO
	829
	830	i = i_start-1
	831	DO j = j_start, MIN( jte, jde-1 )
	832	DO k = kts, kte
	833	hat(i,k,j) = w(i,k,j) / msftx(i,j)
	834	END DO
	835	END DO
	836
	837	j = j_start-1
	838	DO k = kts, kte
	839	DO i = i_start, MIN( ite, ide-1 )
	840	hat(i,k,j) = w(i,k,j) / msftx(i,j)
	841	END DO
	842	END DO
	843
	844	! QUESTION: What is this for?
	845
	846	DO j = j_start, j_end
	847	DO k = kts, ktf
	848	DO i = i_start, i_end
	849	hatavg(i,k,j) = 0.25 * ( &
	850	hat(i,k ,j ) + &
	851	hat(i,k+1,j ) + &
	852	hat(i,k ,j-1) + &
	853	hat(i,k+1,j-1) )
	854	END DO
	855	END DO
	856	END DO
	857
	858	! Calculate dw/dy and store in tmp1.
	859
	860	DO j = j_start, j_end
	861	DO k = kts+1, ktf
	862	DO i = i_start, i_end
	863	tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zy(i,k,j) * &
	864	0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
	865	END DO
	866	END DO
	867	END DO
	868
	869	! End calculation of dw/dy.
	870	!-----------------------------------------------------------------------
	871
	872	!-----------------------------------------------------------------------
	873	! Add the first term (dw/dy) to defor23 (dw/dy+dv/dz) at vorticity
	874	! points.
	875
	876	DO j = j_start, j_end
	877	DO k = kts+1, ktf
	878	DO i = i_start, i_end
	879	defor23(i,k,j) = mm(i,j) * ( &
	880	rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) )
	881	END DO
	882	END DO
	883	END DO
	884
	885	DO j = j_start, j_end
	886	DO i = i_start, i_end
	887	defor23(i,kts,j ) = 0.0
	888	defor23(i,ktf+1,j) = 0.0
	889	END DO
	890	END DO
	891
	892	! End addition of the first term to defor23.
	893	!-----------------------------------------------------------------------
	894
	895	!-----------------------------------------------------------------------
	896	! Calculate dv/dz.
	897
	898	IF ( config_flags%mix_full_fields ) THEN
	899
	900	DO j = j_start, j_end
	901	DO k = kts+1, ktf
	902	DO i = i_start, i_end
	903	tmp1(i,k,j) = ( v(i,k,j) - v(i,k-1,j) ) * &
	904	0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
	905	END DO
	906	END DO
	907	END DO
	908
	909	ELSE
	910
	911	DO j = j_start, j_end
	912	DO k = kts+1, ktf
	913	DO i = i_start, i_end
	914	tmp1(i,k,j) = ( v(i,k,j) - v_base(k) - v(i,k-1,j) + v_base(k-1) ) * &
	915	0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
	916	END DO
	917	END DO
	918	END DO
	919
	920	END IF
	921
	922	! End calculation of dv/dz.
	923	!-----------------------------------------------------------------------
	924
	925	!-----------------------------------------------------------------------
	926	! Add the second term (dv/dz) to defor23 (dw/dy+dv/dz) at vorticity
	927	! points.
	928
	929	! Add tmp1 to defor23.
	930
	931	DO j = j_start, j_end
	932	DO k = kts+1, ktf
	933	DO i = i_start, i_end
	934	defor23(i,k,j) = defor23(i,k,j) + tmp1(i,k,j)
	935	END DO
	936	END DO
	937	END DO
	938
	939	! End addition of the second term to defor23.
	940	!-----------------------------------------------------------------------
	941
	942	!-----------------------------------------------------------------------
	943	! Update the boundary for defor13 and defor23 (might need to change
	944	! later).
	945
	946	IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1) THEN
	947	DO j = jts, jte
	948	DO k = kts, kte
	949	defor13(ids,k,j) = defor13(ids+1,k,j)
	950	defor23(ids,k,j) = defor23(ids+1,k,j)
	951	END DO
	952	END DO
	953	END IF
	954
	955	IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
	956	DO k = kts, kte
	957	DO i = its, ite
	958	defor13(i,k,jds) = defor13(i,k,jds+1)
	959	defor23(i,k,jds) = defor23(i,k,jds+1)
	960	END DO
	961	END DO
	962	END IF
	963
	964	IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
	965	DO j = jts, jte
	966	DO k = kts, kte
	967	defor13(ide,k,j) = defor13(ide-1,k,j)
	968	defor23(ide,k,j) = defor23(ide-1,k,j)
	969	END DO
	970	END DO
	971	END IF
	972
	973	IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
	974	DO k = kts, kte
	975	DO i = its, ite
	976	defor13(i,k,jde) = defor13(i,k,jde-1)
	977	defor23(i,k,jde) = defor23(i,k,jde-1)
	978	END DO
	979	END DO
	980	END IF
	981
	982	! End update of boundary for defor13 and defor23.
	983	!-----------------------------------------------------------------------
	984
	985	! The second three (defor12, defor13, defor23) of six deformation terms
	986	! are now calculated at vorticity points.
	987	!=======================================================================
	988
	989	END SUBROUTINE cal_deform_and_div
	990
	991	!=======================================================================
	992	!=======================================================================
	993
	994	SUBROUTINE calculate_km_kh( config_flags, dt, &
	995	dampcoef, zdamp, damp_opt, &
	996	xkmh, xkmv, xkhh, xkhv, &
	997	BN2, khdif, kvdif, div, &
	998	defor11, defor22, defor33, &
	999	defor12, defor13, defor23, &
	1000	tke, p8w, t8w, theta, t, p, moist, &
	1001	dn, dnw, dx, dy, rdz, rdzw, isotropic, &
	1002	n_moist, cf1, cf2, cf3, warm_rain, &
	1003	mix_upper_bound, &
	1004	msftx, msfty, &
	1005	ids, ide, jds, jde, kds, kde, &
	1006	ims, ime, jms, jme, kms, kme, &
	1007	its, ite, jts, jte, kts, kte )
	1008
	1009	! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
	1010	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	1011	! ... ...
	1012
	1013	! Purpose: This routine calculates exchange coefficients for the TKE
	1014	! scheme.
	1015
	1016	! References: Klemp and Wilhelmson (JAS 1978)
	1017	! Deardorff (B-L Meteor 1980)
	1018	! Chen and Dudhia (NCAR WRF physics report 2000)
	1019
	1020	!-----------------------------------------------------------------------
	1021	! Begin declarations.
	1022
	1023	IMPLICIT NONE
	1024
	1025	TYPE( grid_config_rec_type ), INTENT( IN ) &
	1026	:: config_flags
	1027
	1028	INTEGER, INTENT( IN ) &
	1029	:: n_moist, damp_opt, isotropic, &
	1030	ids, ide, jds, jde, kds, kde, &
	1031	ims, ime, jms, jme, kms, kme, &
	1032	its, ite, jts, jte, kts, kte
	1033
	1034	LOGICAL, INTENT( IN ) &
	1035	:: warm_rain
	1036
	1037	REAL, INTENT( IN ) &
	1038	:: dx, dy, zdamp, dt, dampcoef, cf1, cf2, cf3, khdif, kvdif
	1039
	1040	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
	1041	:: dnw, dn
	1042
	1043	REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist ), INTENT( INOUT ) &
	1044	:: moist
	1045
	1046	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	1047	:: xkmv, xkmh, xkhv, xkhh, BN2
	1048
	1049	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT( IN ) &
	1050	:: defor11, defor22, defor33, defor12, defor13, defor23, &
	1051	div, rdz, rdzw, p8w, t8w, theta, t, p
	1052
	1053	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	1054	:: tke
	1055
	1056	REAL, INTENT( IN ) &
	1057	:: mix_upper_bound
	1058
	1059	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	1060	:: msftx, msfty
	1061
	1062	! Local variables.
	1063
	1064	INTEGER &
	1065	:: i_start, i_end, j_start, j_end, ktf, i, j, k
	1066
	1067	! End declarations.
	1068	!-----------------------------------------------------------------------
	1069
	1070	ktf = MIN( kte, kde-1 )
	1071	i_start = its
	1072	i_end = MIN( ite, ide-1 )
	1073	j_start = jts
	1074	j_end = MIN( jte, jde-1 )
	1075
	1076	CALL calculate_N2( config_flags, BN2, moist, &
	1077	theta, t, p, p8w, t8w, &
	1078	dnw, dn, rdz, rdzw, &
	1079	n_moist, cf1, cf2, cf3, warm_rain, &
	1080	ids, ide, jds, jde, kds, kde, &
	1081	ims, ime, jms, jme, kms, kme, &
	1082	its, ite, jts, jte, kts, kte )
	1083
	1084	! Select a scheme for calculating diffusion coefficients.
	1085
	1086	km_coef: SELECT CASE( config_flags%km_opt )
	1087
	1088	CASE (1)
	1089	CALL isotropic_km( config_flags, xkmh, xkmv, &
	1090	xkhh, xkhv, khdif, kvdif, &
	1091	ids, ide, jds, jde, kds, kde, &
	1092	ims, ime, jms, jme, kms, kme, &
	1093	its, ite, jts, jte, kts, kte )
	1094	CASE (2)
	1095	CALL tke_km( config_flags, xkmh, xkmv, &
	1096	xkhh, xkhv, BN2, tke, p8w, t8w, theta, &
	1097	rdz, rdzw, dx, dy, dt, isotropic, &
	1098	mix_upper_bound, msftx, msfty, &
	1099	ids, ide, jds, jde, kds, kde, &
	1100	ims, ime, jms, jme, kms, kme, &
	1101	its, ite, jts, jte, kts, kte )
	1102	CASE (3)
	1103	CALL smag_km( config_flags, xkmh, xkmv, &
	1104	xkhh, xkhv, BN2, div, &
	1105	defor11, defor22, defor33, &
	1106	defor12, defor13, defor23, &
	1107	rdzw, dx, dy, dt, isotropic, &
	1108	mix_upper_bound, msftx, msfty, &
	1109	ids, ide, jds, jde, kds, kde, &
	1110	ims, ime, jms, jme, kms, kme, &
	1111	its, ite, jts, jte, kts, kte )
	1112	CASE (4)
	1113	CALL smag2d_km( config_flags, xkmh, xkmv, &
	1114	xkhh, xkhv, defor11, defor22, defor12, &
	1115	rdzw, dx, dy, msftx, msfty, &
	1116	ids, ide, jds, jde, kds, kde, &
	1117	ims, ime, jms, jme, kms, kme, &
	1118	its, ite, jts, jte, kts, kte )
	1119	CASE DEFAULT
	1120	CALL wrf_error_fatal( 'Please choose diffusion coefficient scheme' )
	1121
	1122	END SELECT km_coef
	1123
	1124	IF ( damp_opt .eq. 1 ) THEN
	1125	CALL cal_dampkm( config_flags, xkmh, xkhh, xkmv, xkhv, &
	1126	dx, dy, dt, dampcoef, rdz, rdzw, zdamp, &
	1127	msftx, msfty, &
	1128	ids, ide, jds, jde, kds, kde, &
	1129	ims, ime, jms, jme, kms, kme, &
	1130	its, ite, jts, jte, kts, kte )
	1131	END IF
	1132
	1133	END SUBROUTINE calculate_km_kh
	1134
	1135	!=======================================================================
	1136
	1137	SUBROUTINE cal_dampkm( config_flags,xkmh,xkhh,xkmv,xkhv, &
	1138	dx,dy,dt,dampcoef, &
	1139	rdz, rdzw ,zdamp, &
	1140	msftx, msfty, &
	1141	ids,ide, jds,jde, kds,kde, &
	1142	ims,ime, jms,jme, kms,kme, &
	1143	its,ite, jts,jte, kts,kte )
	1144
	1145	!-----------------------------------------------------------------------
	1146	! Begin declarations.
	1147
	1148	IMPLICIT NONE
	1149
	1150	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
	1151
	1152	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	1153	ims, ime, jms, jme, kms, kme, &
	1154	its, ite, jts, jte, kts, kte
	1155
	1156	REAL , INTENT(IN ) :: zdamp,dx,dy,dt,dampcoef
	1157
	1158
	1159	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: xkmh , &
	1160	xkhh , &
	1161	xkmv , &
	1162	xkhv
	1163
	1164	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdz, &
	1165	rdzw
	1166
	1167	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
	1168	msfty
	1169	! LOCAL VARS
	1170
	1171	INTEGER :: i_start, i_end, j_start, j_end, ktf, ktfm1, i, j, k
	1172	REAL :: kmmax,kmmvmax,degrad90,dz,tmp
	1173	REAL :: ds
	1174	REAL , DIMENSION( its:ite ) :: deltaz
	1175	REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: dampk,dampkv
	1176
	1177	! End declarations.
	1178	!-----------------------------------------------------------------------
	1179
	1180	ktf = min(kte,kde-1)
	1181	ktfm1 = ktf-1
	1182
	1183	i_start = its
	1184	i_end = MIN(ite,ide-1)
	1185	j_start = jts
	1186	j_end = MIN(jte,jde-1)
	1187
	1188	! keep upper damping diffusion away from relaxation zones at boundaries if used
	1189	IF(config_flags%specified .OR. config_flags%nested)THEN
	1190	i_start = MAX(i_start,ids+config_flags%spec_bdy_width-1)
	1191	i_end = MIN(i_end,ide-config_flags%spec_bdy_width)
	1192	j_start = MAX(j_start,jds+config_flags%spec_bdy_width-1)
	1193	j_end = MIN(j_end,jde-config_flags%spec_bdy_width)
	1194	ENDIF
	1195
	1196	kmmax=dx*dx/dt
	1197	degrad90=DEGRAD*90.
	1198	DO j = j_start, j_end
	1199
	1200	k=ktf
	1201	DO i = i_start, i_end
	1202	! Unmodified dx used above may produce very large diffusivities
	1203	! when msftx is very large. And the above formula ignores the fact
	1204	! that dy may now be different from dx as well. Let's fix that by
	1205	! defining a "ds" as the minimum of the "real-space" (physical
	1206	! distance) values of dx and dy, and then using that smallest value
	1207	! to calculate a point-by-point kmmax
	1208	ds = MIN(dx/msftx(i,j),dy/msfty(i,j))
	1209	kmmax=ds*ds/dt
	1210
	1211	! deltaz(i)=0.5*dnw(k)/zeta_z(i,j)
	1212	! dz=dnw(k)/zeta_z(i,j)
	1213	dz = 1./rdzw(i,k,j)
	1214	deltaz(i) = 0.5*dz
	1215
	1216	kmmvmax=dz*dz/dt
	1217	tmp=min(deltaz(i)/zdamp,1.)
	1218	dampk(i,k,j)=cos(degrad90tmp)cos(degrad90tmp)kmmax*dampcoef
	1219	dampkv(i,k,j)=cos(degrad90tmp)cos(degrad90tmp)kmmvmax*dampcoef
	1220	! set upper limit on vertical K (based on horizontal K)
	1221	dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j))
	1222
	1223	ENDDO
	1224
	1225	DO k = ktfm1,kts,-1
	1226	DO i = i_start, i_end
	1227	! Unmodified dx used above may produce very large diffusivities
	1228	! when msftx is very large. And the above formula ignores the fact
	1229	! that dy may now be different from dx as well. Let's fix that by
	1230	! defining a "ds" as the minimum of the "real-space" (physical
	1231	! distance) values of dx and dy, and then using that smallest value
	1232	! to calculate a point-by-point kmmax
	1233	ds = MIN(dx/msftx(i,j),dy/msfty(i,j))
	1234	kmmax=ds*ds/dt
	1235
	1236	! deltaz(i)=deltaz(i)+dn(k)/zeta_z(i,j)
	1237	! dz=dnw(k)/zeta_z(i,j)
	1238	dz = 1./rdz(i,k,j)
	1239	deltaz(i) = deltaz(i) + dz
	1240	dz = 1./rdzw(i,k,j)
	1241
	1242	kmmvmax=dz*dz/dt
	1243	tmp=min(deltaz(i)/zdamp,1.)
	1244	dampk(i,k,j)=cos(degrad90tmp)cos(degrad90tmp)kmmax*dampcoef
	1245	dampkv(i,k,j)=cos(degrad90tmp)cos(degrad90tmp)kmmvmax*dampcoef
	1246	! set upper limit on vertical K (based on horizontal K)
	1247	dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j))
	1248	ENDDO
	1249	ENDDO
	1250
	1251	ENDDO
	1252
	1253	DO j = j_start, j_end
	1254	DO k = kts,ktf
	1255	DO i = i_start, i_end
	1256	xkmh(i,k,j)=max(xkmh(i,k,j),dampk(i,k,j))
	1257	xkhh(i,k,j)=max(xkhh(i,k,j),dampk(i,k,j))
	1258	xkmv(i,k,j)=max(xkmv(i,k,j),dampkv(i,k,j))
	1259	xkhv(i,k,j)=max(xkhv(i,k,j),dampkv(i,k,j))
	1260	ENDDO
	1261	ENDDO
	1262	ENDDO
	1263
	1264	END SUBROUTINE cal_dampkm
	1265
	1266	!=======================================================================
	1267	!=======================================================================
	1268
	1269	SUBROUTINE calculate_N2( config_flags, BN2, moist, &
	1270	theta, t, p, p8w, t8w, &
	1271	dnw, dn, rdz, rdzw, &
	1272	n_moist, cf1, cf2, cf3, warm_rain, &
	1273	ids, ide, jds, jde, kds, kde, &
	1274	ims, ime, jms, jme, kms, kme, &
	1275	its, ite, jts, jte, kts, kte )
	1276
	1277	!-----------------------------------------------------------------------
	1278	! Begin declarations.
	1279
	1280	IMPLICIT NONE
	1281
	1282	TYPE( grid_config_rec_type ), INTENT( IN ) &
	1283	:: config_flags
	1284
	1285	INTEGER, INTENT( IN ) &
	1286	:: n_moist, &
	1287	ids, ide, jds, jde, kds, kde, &
	1288	ims, ime, jms, jme, kms, kme, &
	1289	its, ite, jts, jte, kts, kte
	1290
	1291	LOGICAL, INTENT( IN ) &
	1292	:: warm_rain
	1293
	1294	REAL, INTENT( IN ) &
	1295	:: cf1, cf2, cf3
	1296
	1297	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	1298	:: BN2
	1299
	1300	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	1301	:: rdz, rdzw, theta, t, p, p8w, t8w
	1302
	1303	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
	1304	:: dnw, dn
	1305
	1306	REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), INTENT( INOUT ) &
	1307	:: moist
	1308
	1309	! Local variables.
	1310
	1311	INTEGER &
	1312	:: i, j, k, ktf, ispe, ktes1, ktes2, &
	1313	i_start, i_end, j_start, j_end
	1314
	1315	REAL &
	1316	:: coefa, thetaep1, thetaem1, qc_cr, es, tc, qlpqi, qsw, qsi, &
	1317	tmpdz, xlvqv, thetaesfc, thetasfc, qvtop, qvsfc, thetatop, thetaetop
	1318
	1319	REAL, DIMENSION( its:ite, jts:jte ) &
	1320	:: tmp1sfc, tmp1top
	1321
	1322	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
	1323	:: tmp1, qvs, qctmp
	1324
	1325	! End declarations.
	1326	!-----------------------------------------------------------------------
	1327
	1328	qc_cr = 0.00001 ! in Kg/Kg
	1329
	1330	ktf = MIN( kte, kde-1 )
	1331	ktes1 = kte-1
	1332	ktes2 = kte-2
	1333
	1334	i_start = its
	1335	i_end = MIN( ite, ide-1 )
	1336	j_start = jts
	1337	j_end = MIN( jte, jde-1 )
	1338
	1339	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	1340	config_flags%nested) i_start = MAX( ids+1, its )
	1341	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	1342	config_flags%nested) i_end = MIN( ide-2, ite )
	1343	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	1344	config_flags%nested) j_start = MAX( jds+1, jts )
	1345	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	1346	config_flags%nested) j_end = MIN( jde-2 ,jte )
	1347	IF ( config_flags%periodic_x ) i_start = its
	1348	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	1349
	1350	IF ( P_QC .GT. PARAM_FIRST_SCALAR) THEN
	1351	DO j = j_start, j_end
	1352	DO k = kts, ktf
	1353	DO i = i_start, i_end
	1354	qctmp(i,k,j) = moist(i,k,j,P_QC)
	1355	END DO
	1356	END DO
	1357	END DO
	1358	ELSE
	1359	DO j = j_start, j_end
	1360	DO k = kts, ktf
	1361	DO i = i_start, i_end
	1362	qctmp(i,k,j) = 0.0
	1363	END DO
	1364	END DO
	1365	END DO
	1366	END IF
	1367
	1368	DO j = jts, jte
	1369	DO k = kts, kte
	1370	DO i = its, ite
	1371	tmp1(i,k,j) = 0.0
	1372	END DO
	1373	END DO
	1374	END DO
	1375
	1376	DO j = jts,jte
	1377	DO i = its,ite
	1378	tmp1sfc(i,j) = 0.0
	1379	tmp1top(i,j) = 0.0
	1380	END DO
	1381	END DO
	1382
	1383	DO ispe = PARAM_FIRST_SCALAR, n_moist
	1384	IF ( ispe .EQ. P_QV .OR. ispe .EQ. P_QC .OR. ispe .EQ. P_QI) THEN
	1385	DO j = j_start, j_end
	1386	DO k = kts, ktf
	1387	DO i = i_start, i_end
	1388	tmp1(i,k,j) = tmp1(i,k,j) + moist(i,k,j,ispe)
	1389	END DO
	1390	END DO
	1391	END DO
	1392
	1393	DO j = j_start, j_end
	1394	DO i = i_start, i_end
	1395	tmp1sfc(i,j) = tmp1sfc(i,j) + &
	1396	cf1 * moist(i,1,j,ispe) + &
	1397	cf2 * moist(i,2,j,ispe) + &
	1398	cf3 * moist(i,3,j,ispe)
	1399	tmp1top(i,j) = tmp1top(i,j) + &
	1400	moist(i,ktes1,j,ispe) + &
	1401	( moist(i,ktes1,j,ispe) - moist(i,ktes2,j,ispe) ) * &
	1402	0.5 * dnw(ktes1) / dn(ktes1)
	1403	END DO
	1404	END DO
	1405	END IF
	1406	END DO
	1407
	1408	! Calculate saturation mixing ratio.
	1409
	1410	DO j = j_start, j_end
	1411	DO k = kts, ktf
	1412	DO i = i_start, i_end
	1413	tc = t(i,k,j) - SVPT0
	1414	es = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t(i,k,j) - SVP3 ) )
	1415	qvs(i,k,j) = EP_2 * es / ( p(i,k,j) - es )
	1416	END DO
	1417	END DO
	1418	END DO
	1419
	1420	DO j = j_start, j_end
	1421	DO k = kts+1, ktf-1
	1422	DO i = i_start, i_end
	1423	tmpdz = 1.0 / rdz(i,k,j) + 1.0 / rdz(i,k+1,j)
	1424	IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN
	1425	xlvqv = XLV * moist(i,k,j,P_QV)
	1426	coefa = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / &
	1427	( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) / &
	1428	theta(i,k,j)
	1429	thetaep1 = theta(i,k+1,j) * &
	1430	( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) )
	1431	thetaem1 = theta(i,k-1,j) * &
	1432	( 1.0 + XLV * qvs(i,k-1,j) / Cp / t(i,k-1,j) )
	1433	BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaem1 ) / tmpdz - &
	1434	( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz )
	1435	ELSE
	1436	BN2(i,k,j) = g * ( (theta(i,k+1,j) - theta(i,k-1,j) ) / &
	1437	theta(i,k,j) / tmpdz + &
	1438	1.61 * ( moist(i,k+1,j,P_QV) - moist(i,k-1,j,P_QV) ) / &
	1439	tmpdz - &
	1440	( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz )
	1441	ENDIF
	1442	END DO
	1443	END DO
	1444	END DO
	1445
	1446	k = kts
	1447	DO j = j_start, j_end
	1448	DO i = i_start, i_end
	1449	tmpdz = 1.0 / rdz(i,k+1,j) + 0.5 / rdzw(i,k,j)
	1450	thetasfc = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp )
	1451	IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN
	1452	qvsfc = cf1 * qvs(i,1,j) + &
	1453	cf2 * qvs(i,2,j) + &
	1454	cf3 * qvs(i,3,j)
	1455	xlvqv = XLV * moist(i,k,j,P_QV)
	1456	coefa = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / &
	1457	( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) / &
	1458	theta(i,k,j)
	1459	thetaep1 = theta(i,k+1,j) * &
	1460	( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) )
	1461	thetaesfc = thetasfc * &
	1462	( 1.0 + XLV * qvsfc / Cp / t8w(i,kts,j) )
	1463	BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaesfc ) / tmpdz - &
	1464	( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
	1465	ELSE
	1466	qvsfc = cf1 * moist(i,1,j,P_QV) + &
	1467	cf2 * moist(i,2,j,P_QV) + &
	1468	cf3 * moist(i,3,j,P_QV)
	1469	! BN2(i,k,j) = g * ( ( theta(i,k+1,j) - thetasfc ) / &
	1470	! theta(i,k,j) / tmpdz + &
	1471	! 1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) / &
	1472	! tmpdz - &
	1473	! ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
	1474	!...... MARTA: change in computation of BN2 at the surface, WCS 040331
	1475
	1476	tmpdz= 1./rdzw(i,k,j) ! controlare come calcola rdzw
	1477	BN2(i,k,j) = g * ( ( theta(i,k+1,j) - theta(i,k,j)) / &
	1478	theta(i,k,j) / tmpdz + &
	1479	1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) / &
	1480	tmpdz - &
	1481	( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
	1482	! end of MARTA/WCS change
	1483
	1484	ENDIF
	1485	END DO
	1486	END DO
	1487
	1488
	1489	!...... MARTA: change in computation of BN2 at the top, WCS 040331
	1490	DO j = j_start, j_end
	1491	DO i = i_start, i_end
	1492	BN2(i,ktf,j)=BN2(i,ktf-1,j)
	1493	END DO
	1494	END DO
	1495	! end of MARTA/WCS change
	1496
	1497	END SUBROUTINE calculate_N2
	1498
	1499	!=======================================================================
	1500	!=======================================================================
	1501
	1502	SUBROUTINE isotropic_km( config_flags, &
	1503	xkmh,xkmv,xkhh,xkhv,khdif,kvdif, &
	1504	ids,ide, jds,jde, kds,kde, &
	1505	ims,ime, jms,jme, kms,kme, &
	1506	its,ite, jts,jte, kts,kte )
	1507
	1508	!-----------------------------------------------------------------------
	1509	! Begin declarations.
	1510
	1511	IMPLICIT NONE
	1512
	1513	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
	1514
	1515	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	1516	ims, ime, jms, jme, kms, kme, &
	1517	its, ite, jts, jte, kts, kte
	1518
	1519	REAL , INTENT(IN ) :: khdif,kvdif
	1520
	1521	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
	1522	xkmv, &
	1523	xkhh, &
	1524	xkhv
	1525	! LOCAL VARS
	1526
	1527	INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
	1528	REAL :: khdif3,kvdif3
	1529
	1530	! End declarations.
	1531	!-----------------------------------------------------------------------
	1532
	1533	ktf = kte
	1534
	1535	i_start = its
	1536	i_end = MIN(ite,ide-1)
	1537	j_start = jts
	1538	j_end = MIN(jte,jde-1)
	1539
	1540	! khdif3=khdif*3.
	1541	! kvdif3=kvdif*3.
	1542	khdif3=khdif/prandtl
	1543	kvdif3=kvdif/prandtl
	1544
	1545	DO j = j_start, j_end
	1546	DO k = kts, ktf
	1547	DO i = i_start, i_end
	1548	xkmh(i,k,j)=khdif
	1549	xkmv(i,k,j)=kvdif
	1550	xkhh(i,k,j)=khdif3
	1551	xkhv(i,k,j)=kvdif3
	1552	ENDDO
	1553	ENDDO
	1554	ENDDO
	1555
	1556	END SUBROUTINE isotropic_km
	1557
	1558	!=======================================================================
	1559	!=======================================================================
	1560
	1561	SUBROUTINE smag_km( config_flags,xkmh,xkmv,xkhh,xkhv,BN2, &
	1562	div,defor11,defor22,defor33,defor12, &
	1563	defor13,defor23, &
	1564	rdzw,dx,dy,dt,isotropic, &
	1565	mix_upper_bound, msftx, msfty, &
	1566	ids,ide, jds,jde, kds,kde, &
	1567	ims,ime, jms,jme, kms,kme, &
	1568	its,ite, jts,jte, kts,kte )
	1569
	1570	!-----------------------------------------------------------------------
	1571	! Begin declarations.
	1572
	1573	IMPLICIT NONE
	1574
	1575	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
	1576
	1577	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	1578	ims, ime, jms, jme, kms, kme, &
	1579	its, ite, jts, jte, kts, kte
	1580
	1581	INTEGER , INTENT(IN ) :: isotropic
	1582	REAL , INTENT(IN ) :: dx, dy, dt, mix_upper_bound
	1583
	1584
	1585	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: BN2, &
	1586	rdzw
	1587
	1588	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
	1589	xkmv, &
	1590	xkhh, &
	1591	xkhv
	1592
	1593	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
	1594	defor11, &
	1595	defor22, &
	1596	defor33, &
	1597	defor12, &
	1598	defor13, &
	1599	defor23, &
	1600	div
	1601	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
	1602	msfty
	1603	! LOCAL VARS
	1604
	1605	INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
	1606	REAL :: deltas, tmp, pr, mlen_h, mlen_v, c_s
	1607
	1608	REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2
	1609
	1610	! End declarations.
	1611	!-----------------------------------------------------------------------
	1612
	1613	ktf = min(kte,kde-1)
	1614
	1615	i_start = its
	1616	i_end = MIN(ite,ide-1)
	1617	j_start = jts
	1618	j_end = MIN(jte,jde-1)
	1619
	1620	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	1621	config_flags%nested) i_start = MAX(ids+1,its)
	1622	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	1623	config_flags%nested) i_end = MIN(ide-2,ite)
	1624	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	1625	config_flags%nested) j_start = MAX(jds+1,jts)
	1626	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	1627	config_flags%nested) j_end = MIN(jde-2,jte)
	1628	IF ( config_flags%periodic_x ) i_start = its
	1629	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	1630
	1631	pr = prandtl
	1632	c_s = config_flags%c_s
	1633
	1634	do j=j_start,j_end
	1635	do k=kts,ktf
	1636	do i=i_start,i_end
	1637	def2(i,k,j)=0.5(defor11(i,k,j)defor11(i,k,j) + &
	1638	defor22(i,k,j)*defor22(i,k,j) + &
	1639	defor33(i,k,j)*defor33(i,k,j))
	1640	enddo
	1641	enddo
	1642	enddo
	1643
	1644	do j=j_start,j_end
	1645	do k=kts,ktf
	1646	do i=i_start,i_end
	1647	tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ &
	1648	defor12(i+1,k,j)+defor12(i+1,k,j+1))
	1649	def2(i,k,j)=def2(i,k,j)+tmp*tmp
	1650	enddo
	1651	enddo
	1652	enddo
	1653
	1654	do j=j_start,j_end
	1655	do k=kts,ktf
	1656	do i=i_start,i_end
	1657	tmp=0.25*(defor13(i ,k+1,j)+defor13(i ,k,j)+ &
	1658	defor13(i+1,k+1,j)+defor13(i+1,k,j))
	1659	def2(i,k,j)=def2(i,k,j)+tmp*tmp
	1660	enddo
	1661	enddo
	1662	enddo
	1663
	1664	do j=j_start,j_end
	1665	do k=kts,ktf
	1666	do i=i_start,i_end
	1667	tmp=0.25*(defor23(i,k+1,j )+defor23(i,k,j )+ &
	1668	defor23(i,k+1,j+1)+defor23(i,k,j+1))
	1669	def2(i,k,j)=def2(i,k,j)+tmp*tmp
	1670	enddo
	1671	enddo
	1672	enddo
	1673	!
	1674	IF (isotropic .EQ. 0) THEN
	1675	DO j = j_start, j_end
	1676	DO k = kts, ktf
	1677	DO i = i_start, i_end
	1678	mlen_h=sqrt(dx/msftx(i,j) * dy/msfty(i,j))
	1679	mlen_v= 1./rdzw(i,k,j)
	1680	tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
	1681	tmp=tmp**0.5
	1682	xkmh(i,k,j)=max(c_sc_smlen_hmlen_htmp, 1.0E-6mlen_hmlen_h )
	1683	xkmh(i,k,j)=min(xkmh(i,k,j), mix_upper_bound * mlen_h * mlen_h / dt )
	1684	xkmv(i,k,j)=max(c_sc_smlen_vmlen_vtmp, 1.0E-6mlen_vmlen_v )
	1685	xkmv(i,k,j)=min(xkmv(i,k,j), mix_upper_bound * mlen_v * mlen_v / dt )
	1686	xkhh(i,k,j)=xkmh(i,k,j)/pr
	1687	xkhh(i,k,j)=min(xkhh(i,k,j), mix_upper_bound * mlen_h * mlen_h / dt )
	1688	xkhv(i,k,j)=xkmv(i,k,j)/pr
	1689	xkhv(i,k,j)=min(xkhv(i,k,j), mix_upper_bound * mlen_v * mlen_v / dt )
	1690	ENDDO
	1691	ENDDO
	1692	ENDDO
	1693	ELSE
	1694	DO j = j_start, j_end
	1695	DO k = kts, ktf
	1696	DO i = i_start, i_end
	1697	deltas=(dx/msftx(i,j) * dy/msfty(i,j)/rdzw(i,k,j))**0.33333333
	1698	tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
	1699	tmp=tmp**0.5
	1700	xkmh(i,k,j)=max(c_sc_sdeltasdeltastmp, 1.0E-6deltasdeltas )
	1701	xkmh(i,k,j)=min(xkmh(i,k,j), mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt )
	1702	xkmv(i,k,j)=xkmh(i,k,j)
	1703	xkmv(i,k,j)=min(xkmv(i,k,j), mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt )
	1704	xkhh(i,k,j)=xkmh(i,k,j)/pr
	1705	xkhh(i,k,j)=min(xkhh(i,k,j), mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt )
	1706	xkhv(i,k,j)=xkmv(i,k,j)/pr
	1707	xkhv(i,k,j)=min(xkhv(i,k,j), mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt )
	1708	ENDDO
	1709	ENDDO
	1710	ENDDO
	1711	ENDIF
	1712
	1713	END SUBROUTINE smag_km
	1714
	1715	!=======================================================================
	1716	!=======================================================================
	1717
	1718	SUBROUTINE smag2d_km( config_flags,xkmh,xkmv,xkhh,xkhv, &
	1719	defor11,defor22,defor12, &
	1720	rdzw,dx,dy,msftx, msfty, &
	1721	ids,ide, jds,jde, kds,kde, &
	1722	ims,ime, jms,jme, kms,kme, &
	1723	its,ite, jts,jte, kts,kte )
	1724
	1725	!-----------------------------------------------------------------------
	1726	! Begin declarations.
	1727
	1728	IMPLICIT NONE
	1729
	1730	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
	1731
	1732	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	1733	ims, ime, jms, jme, kms, kme, &
	1734	its, ite, jts, jte, kts, kte
	1735
	1736	REAL , INTENT(IN ) :: dx, dy
	1737
	1738
	1739	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rdzw
	1740
	1741	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
	1742	xkmv, &
	1743	xkhh, &
	1744	xkhv
	1745
	1746	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
	1747	defor11, &
	1748	defor22, &
	1749	defor12
	1750
	1751	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
	1752	msfty
	1753	! LOCAL VARS
	1754
	1755	INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
	1756	REAL :: deltas, tmp, pr, mlen_h, c_s
	1757
	1758	REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2
	1759
	1760	! End declarations.
	1761	!-----------------------------------------------------------------------
	1762
	1763	ktf = min(kte,kde-1)
	1764
	1765	i_start = its
	1766	i_end = MIN(ite,ide-1)
	1767	j_start = jts
	1768	j_end = MIN(jte,jde-1)
	1769
	1770	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	1771	config_flags%nested) i_start = MAX(ids+1,its)
	1772	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	1773	config_flags%nested) i_end = MIN(ide-2,ite)
	1774	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	1775	config_flags%nested) j_start = MAX(jds+1,jts)
	1776	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	1777	config_flags%nested) j_end = MIN(jde-2,jte)
	1778	IF ( config_flags%periodic_x ) i_start = its
	1779	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	1780
	1781	pr=prandtl
	1782	c_s = config_flags%c_s
	1783
	1784	do j=j_start,j_end
	1785	do k=kts,ktf
	1786	do i=i_start,i_end
	1787	def2(i,k,j)=0.25((defor11(i,k,j)-defor22(i,k,j))(defor11(i,k,j)-defor22(i,k,j)))
	1788	tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ &
	1789	defor12(i+1,k,j)+defor12(i+1,k,j+1))
	1790	def2(i,k,j)=def2(i,k,j)+tmp*tmp
	1791	enddo
	1792	enddo
	1793	enddo
	1794	!
	1795	DO j = j_start, j_end
	1796	DO k = kts, ktf
	1797	DO i = i_start, i_end
	1798	mlen_h=sqrt(dx/msftx(i,j) * dy/msfty(i,j))
	1799	tmp=sqrt(def2(i,k,j))
	1800	! xkmh(i,k,j)=max(c_sc_smlen_hmlen_htmp, 1.0E-6mlen_hmlen_h )
	1801	xkmh(i,k,j)=c_sc_smlen_hmlen_htmp
	1802	xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h )
	1803	xkmv(i,k,j)=0.
	1804	xkhh(i,k,j)=xkmh(i,k,j)/pr
	1805	xkhv(i,k,j)=0.
	1806	ENDDO
	1807	ENDDO
	1808	ENDDO
	1809
	1810	END SUBROUTINE smag2d_km
	1811
	1812	!=======================================================================
	1813	!=======================================================================
	1814
	1815	SUBROUTINE tke_km( config_flags, xkmh, xkmv, xkhh, xkhv, &
	1816	bn2, tke, p8w, t8w, theta, &
	1817	rdz, rdzw, dx,dy, dt, isotropic, &
	1818	mix_upper_bound, msftx, msfty, &
	1819	ids, ide, jds, jde, kds, kde, &
	1820	ims, ime, jms, jme, kms, kme, &
	1821	its, ite, jts, jte, kts, kte )
	1822
	1823	! History: Sep 2003 Changes by Jason Knievel and George Bryan, NCAR
	1824	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	1825	! ... ...
	1826
	1827	! Purpose: This routine calculates the exchange coefficients for the
	1828	! TKE turbulence parameterization.
	1829
	1830	! References: Klemp and Wilhelmson (JAS 1978)
	1831	! Chen and Dudhia (NCAR WRF physics report 2000)
	1832
	1833	!-----------------------------------------------------------------------
	1834	! Begin declarations.
	1835
	1836	IMPLICIT NONE
	1837
	1838	TYPE( grid_config_rec_type ), INTENT( IN ) &
	1839	:: config_flags
	1840
	1841	INTEGER, INTENT( IN ) &
	1842	:: ids, ide, jds, jde, kds, kde, &
	1843	ims, ime, jms, jme, kms, kme, &
	1844	its, ite, jts, jte, kts, kte
	1845
	1846	INTEGER, INTENT( IN ) :: isotropic
	1847	REAL, INTENT( IN ) &
	1848	:: dx, dy, dt
	1849
	1850	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	1851	:: tke, p8w, t8w, theta, rdz, rdzw, bn2
	1852
	1853	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	1854	:: xkmh, xkmv, xkhh, xkhv
	1855
	1856	REAL, INTENT( IN ) &
	1857	:: mix_upper_bound
	1858
	1859	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
	1860	msfty
	1861	! Local variables.
	1862
	1863	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
	1864	:: l_scale
	1865
	1866	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
	1867	:: dthrdn
	1868
	1869	REAL &
	1870	:: deltas, tmp, mlen_s, mlen_h, mlen_v, tmpdz, &
	1871	thetasfc, thetatop, minkx, pr_inv, pr_inv_h, pr_inv_v, c_k
	1872
	1873	INTEGER &
	1874	:: i_start, i_end, j_start, j_end, ktf, i, j, k
	1875
	1876	REAL, PARAMETER :: tke_seed_value = 1.e-06
	1877	REAL :: tke_seed
	1878	REAL, PARAMETER :: epsilon = 1.e-10
	1879
	1880	! End declarations.
	1881	!-----------------------------------------------------------------------
	1882
	1883	ktf = MIN( kte, kde-1 )
	1884	i_start = its
	1885	i_end = MIN( ite, ide-1 )
	1886	j_start = jts
	1887	j_end = MIN( jte, jde-1 )
	1888
	1889	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	1890	config_flags%nested) i_start = MAX( ids+1, its )
	1891	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	1892	config_flags%nested) i_end = MIN( ide-2, ite )
	1893	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	1894	config_flags%nested) j_start = MAX( jds+1, jts )
	1895	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	1896	config_flags%nested) j_end = MIN( jde-2, jte)
	1897	IF ( config_flags%periodic_x ) i_start = its
	1898	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	1899
	1900	! in the absence of surface drag or a surface heat flux, there
	1901	! is no way to generate tke without pre-existing tke. Use
	1902	! tke_seed if the drag and flux are off.
	1903
	1904	c_k = config_flags%c_k
	1905	tke_seed = tke_seed_value
	1906	if( (config_flags%tke_drag_coefficient .gt. epsilon) .or. &
	1907	(config_flags%tke_heat_flux .gt. epsilon) ) tke_seed = 0.
	1908
	1909	DO j = j_start, j_end
	1910	DO k = kts+1, ktf-1
	1911	DO i = i_start, i_end
	1912	tmpdz = 1.0 / ( rdz(i,k+1,j) + rdz(i,k,j) )
	1913	dthrdn(i,k,j) = ( theta(i,k+1,j) - theta(i,k-1,j) ) / tmpdz
	1914	END DO
	1915	END DO
	1916	END DO
	1917
	1918	k = kts
	1919	DO j = j_start, j_end
	1920	DO i = i_start, i_end
	1921	tmpdz = 1.0 / ( rdzw(i,k+1,j) + rdzw(i,k,j) )
	1922	thetasfc = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp )
	1923	dthrdn(i,k,j) = ( theta(i,k+1,j) - thetasfc ) / tmpdz
	1924	END DO
	1925	END DO
	1926
	1927	k = ktf
	1928	DO j = j_start, j_end
	1929	DO i = i_start, i_end
	1930	tmpdz = 1.0 / rdz(i,k,j) + 0.5 / rdzw(i,k,j)
	1931	thetatop = T8w(i,kde,j) / ( p8w(i,kde,j) / p1000mb )**( R_d / Cp )
	1932	dthrdn(i,k,j) = ( thetatop - theta(i,k-1,j) ) / tmpdz
	1933	END DO
	1934	END DO
	1935
	1936	IF ( isotropic .EQ. 0 ) THEN
	1937	DO j = j_start, j_end
	1938	DO k = kts, ktf
	1939	DO i = i_start, i_end
	1940	mlen_h = SQRT( dx/msftx(i,j) * dy/msfty(i,j) )
	1941	tmp = SQRT( MAX( tke(i,k,j), tke_seed ) )
	1942	deltas = 1.0 / rdzw(i,k,j)
	1943	mlen_v = deltas
	1944	IF ( dthrdn(i,k,j) .GT. 0.) THEN
	1945	mlen_s = 0.76 * tmp / ( ABS( g / theta(i,k,j) * dthrdn(i,k,j) ) )**0.5
	1946	mlen_v = MIN( mlen_v, mlen_s )
	1947	END IF
	1948	xkmh(i,k,j) = MAX( c_k * tmp * mlen_h, 1.0E-6 * mlen_h * mlen_h )
	1949	xkmh(i,k,j) = MIN( xkmh(i,k,j), mix_upper_bound * mlen_h *mlen_h / dt )
	1950	xkmv(i,k,j) = MAX( c_k * tmp * mlen_v, 1.0E-6 * deltas * deltas )
	1951	xkmv(i,k,j) = MIN( xkmv(i,k,j), mix_upper_bound * deltas *deltas / dt )
	1952	pr_inv_h = 1./prandtl
	1953	pr_inv_v = 1.0 + 2.0 * mlen_v / deltas
	1954	xkhh(i,k,j) = xkmh(i,k,j) * pr_inv_h
	1955	xkhv(i,k,j) = xkmv(i,k,j) * pr_inv_v
	1956	END DO
	1957	END DO
	1958	END DO
	1959	ELSE
	1960	CALL calc_l_scale( config_flags, tke, BN2, l_scale, &
	1961	i_start, i_end, ktf, j_start, j_end, &
	1962	dx, dy, rdzw, msftx, msfty, &
	1963	ids, ide, jds, jde, kds, kde, &
	1964	ims, ime, jms, jme, kms, kme, &
	1965	its, ite, jts, jte, kts, kte )
	1966	DO j = j_start, j_end
	1967	DO k = kts, ktf
	1968	DO i = i_start, i_end
	1969	tmp = SQRT( MAX( tke(i,k,j), tke_seed ) )
	1970	deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
	1971	xkmh(i,k,j) = c_k * tmp * l_scale(i,k,j)
	1972	xkmh(i,k,j) = MIN( mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt, xkmh(i,k,j) )
	1973	xkmv(i,k,j) = c_k * tmp * l_scale(i,k,j)
	1974	xkmv(i,k,j) = MIN( mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt , xkmv(i,k,j) )
	1975	pr_inv = 1.0 + 2.0 * l_scale(i,k,j) / deltas
	1976	xkhh(i,k,j) = MIN( mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt, xkmh(i,k,j) * pr_inv )
	1977	xkhv(i,k,j) = MIN( mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt, xkmv(i,k,j) * pr_inv )
	1978	END DO
	1979	END DO
	1980	END DO
	1981	END IF
	1982
	1983	END SUBROUTINE tke_km
	1984
	1985	!=======================================================================
	1986	!=======================================================================
	1987
	1988	SUBROUTINE calc_l_scale( config_flags, tke, BN2, l_scale, &
	1989	i_start, i_end, ktf, j_start, j_end, &
	1990	dx, dy, rdzw, msftx, msfty, &
	1991	ids, ide, jds, jde, kds, kde, &
	1992	ims, ime, jms, jme, kms, kme, &
	1993	its, ite, jts, jte, kts, kte )
	1994
	1995	! History: Sep 2003 Written by Bryan and Knievel, NCAR
	1996
	1997	! Purpose: This routine calculates the length scale, based on stability,
	1998	! for TKE parameterization of subgrid-scale turbulence.
	1999
	2000	!-----------------------------------------------------------------------
	2001	! Begin declarations.
	2002
	2003	IMPLICIT NONE
	2004
	2005	TYPE( grid_config_rec_type ), INTENT( IN ) &
	2006	:: config_flags
	2007
	2008	INTEGER, INTENT( IN ) &
	2009	:: i_start, i_end, ktf, j_start, j_end, &
	2010	ids, ide, jds, jde, kds, kde, &
	2011	ims, ime, jms, jme, kms, kme, &
	2012	its, ite, jts, jte, kts, kte
	2013
	2014	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	2015	:: BN2, tke, rdzw
	2016
	2017	REAL, INTENT( IN ) &
	2018	:: dx, dy
	2019
	2020	REAL, DIMENSION( its:ite, kts:kte, jts:jte ), INTENT( OUT ) &
	2021	:: l_scale
	2022
	2023	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
	2024	msfty
	2025	! Local variables.
	2026
	2027	INTEGER &
	2028	:: i, j, k
	2029
	2030	REAL &
	2031	:: deltas, tmp
	2032
	2033	! End declarations.
	2034	!-----------------------------------------------------------------------
	2035
	2036	DO j = j_start, j_end
	2037	DO k = kts, ktf
	2038	DO i = i_start, i_end
	2039	deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
	2040	l_scale(i,k,j) = deltas
	2041
	2042	IF ( BN2(i,k,j) .gt. 1.0e-6 ) THEN
	2043	tmp = SQRT( MAX( tke(i,k,j), 1.0e-6 ) )
	2044	l_scale(i,k,j) = 0.76 * tmp / SQRT( BN2(i,k,j) )
	2045	l_scale(i,k,j) = MIN( l_scale(i,k,j), deltas)
	2046	l_scale(i,k,j) = MAX( l_scale(i,k,j), 0.001 * deltas )
	2047	END IF
	2048
	2049	END DO
	2050	END DO
	2051	END DO
	2052
	2053	END SUBROUTINE calc_l_scale
	2054
	2055	!=======================================================================
	2056	!=======================================================================
	2057
	2058	SUBROUTINE horizontal_diffusion_2 ( rt_tendf, ru_tendf, rv_tendf, rw_tendf, &
	2059	tke_tendf, &
	2060	moist_tendf, n_moist, &
	2061	chem_tendf, n_chem, &
	2062	scalar_tendf, n_scalar, &
	2063	thp, theta, mu, tke, config_flags, &
	2064	defor11, defor22, defor12, &
	2065	defor13, defor23, div, &
	2066	moist, chem, scalar, &
	2067	msfux, msfuy, msfvx, msfvy, &
	2068	msftx, msfty, xkmh, xkhh,km_opt, &
	2069	rdx, rdy, rdz, rdzw, fnm, fnp, &
	2070	cf1, cf2, cf3, zx, zy, dn, dnw, &
	2071	ids, ide, jds, jde, kds, kde, &
	2072	ims, ime, jms, jme, kms, kme, &
	2073	its, ite, jts, jte, kts, kte )
	2074
	2075	!-----------------------------------------------------------------------
	2076	! Begin declarations.
	2077
	2078	IMPLICIT NONE
	2079
	2080	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	2081
	2082	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	2083	ims, ime, jms, jme, kms, kme, &
	2084	its, ite, jts, jte, kts, kte
	2085
	2086	INTEGER , INTENT(IN ) :: n_moist, n_chem, n_scalar, km_opt
	2087
	2088	REAL , INTENT(IN ) :: cf1, cf2, cf3
	2089
	2090	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	2091	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	2092	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
	2093	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
	2094
	2095	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux, &
	2096	msfuy, &
	2097	msfvx, &
	2098	msfvy, &
	2099	msftx, &
	2100	msfty, &
	2101	mu
	2102
	2103	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::rt_tendf,&
	2104	ru_tendf,&
	2105	rv_tendf,&
	2106	rw_tendf,&
	2107	tke_tendf
	2108
	2109	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
	2110	INTENT(INOUT) :: moist_tendf
	2111
	2112	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
	2113	INTENT(INOUT) :: chem_tendf
	2114
	2115	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar), &
	2116	INTENT(INOUT) :: scalar_tendf
	2117
	2118	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
	2119	INTENT(IN ) :: moist
	2120
	2121	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
	2122	INTENT(IN ) :: chem
	2123
	2124	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
	2125	INTENT(IN ) :: scalar
	2126
	2127	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
	2128	defor22, &
	2129	defor12, &
	2130	defor13, &
	2131	defor23, &
	2132	div, &
	2133	xkmh, &
	2134	xkhh, &
	2135	zx, &
	2136	zy, &
	2137	theta, &
	2138	thp, &
	2139	tke, &
	2140	rdz, &
	2141	rdzw
	2142
	2143
	2144	REAL , INTENT(IN ) :: rdx, &
	2145	rdy
	2146
	2147	! LOCAL VARS
	2148
	2149	INTEGER :: im, ic, is
	2150
	2151	! REAL , DIMENSION(its-1:ite+1, kts:kte, jts-1:jte+1) :: xkhh
	2152
	2153	! End declarations.
	2154	!-----------------------------------------------------------------------
	2155
	2156	!-----------------------------------------------------------------------
	2157	! Call diffusion subroutines.
	2158
	2159	CALL horizontal_diffusion_u_2( ru_tendf, mu, config_flags, &
	2160	defor11, defor12, div, &
	2161	tke(ims,kms,jms), &
	2162	msfux, msfuy, xkmh, rdx, rdy, fnm, fnp, &
	2163	zx, zy, rdzw, &
	2164	ids, ide, jds, jde, kds, kde, &
	2165	ims, ime, jms, jme, kms, kme, &
	2166	its, ite, jts, jte, kts, kte )
	2167
	2168	CALL horizontal_diffusion_v_2( rv_tendf, mu, config_flags, &
	2169	defor12, defor22, div, &
	2170	tke(ims,kms,jms), &
	2171	msfvx, msfvy, xkmh, rdx, rdy, fnm, fnp, &
	2172	zx, zy, rdzw, &
	2173	ids, ide, jds, jde, kds, kde, &
	2174	ims, ime, jms, jme, kms, kme, &
	2175	its, ite, jts, jte, kts, kte )
	2176
	2177	CALL horizontal_diffusion_w_2( rw_tendf, mu, config_flags, &
	2178	defor13, defor23, div, &
	2179	tke(ims,kms,jms), &
	2180	msftx, msfty, xkmh, rdx, rdy, fnm, fnp, &
	2181	zx, zy, rdz, &
	2182	ids, ide, jds, jde, kds, kde, &
	2183	ims, ime, jms, jme, kms, kme, &
	2184	its, ite, jts, jte, kts, kte )
	2185
	2186	CALL horizontal_diffusion_s ( rt_tendf, mu, config_flags, thp, &
	2187	msftx, msfty, msfux, msfuy, &
	2188	msfvx, msfvy, xkhh, rdx, rdy, &
	2189	fnm, fnp, cf1, cf2, cf3, &
	2190	zx, zy, rdz, rdzw, dnw, dn, &
	2191	.false., &
	2192	ids, ide, jds, jde, kds, kde, &
	2193	ims, ime, jms, jme, kms, kme, &
	2194	its, ite, jts, jte, kts, kte )
	2195
	2196	IF (km_opt .eq. 2) &
	2197	CALL horizontal_diffusion_s ( tke_tendf(ims,kms,jms), &
	2198	mu, config_flags, &
	2199	tke(ims,kms,jms), &
	2200	msftx, msfty, msfux, msfuy, &
	2201	msfvx, msfvy, xkhh, rdx, rdy, &
	2202	fnm, fnp, cf1, cf2, cf3, &
	2203	zx, zy, rdz, rdzw, dnw, dn, &
	2204	.true., &
	2205	ids, ide, jds, jde, kds, kde, &
	2206	ims, ime, jms, jme, kms, kme, &
	2207	its, ite, jts, jte, kts, kte )
	2208
	2209	IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
	2210
	2211	moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
	2212
	2213	CALL horizontal_diffusion_s( moist_tendf(ims,kms,jms,im), &
	2214	mu, config_flags, &
	2215	moist(ims,kms,jms,im), &
	2216	msftx, msfty, msfux, msfuy, &
	2217	msfvx, msfvy, xkhh, rdx, rdy, &
	2218	fnm, fnp, cf1, cf2, cf3, &
	2219	zx, zy, rdz, rdzw, dnw, dn, &
	2220	.false., &
	2221	ids, ide, jds, jde, kds, kde, &
	2222	ims, ime, jms, jme, kms, kme, &
	2223	its, ite, jts, jte, kts, kte )
	2224
	2225	ENDDO moist_loop
	2226
	2227	ENDIF
	2228
	2229	IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN
	2230
	2231	chem_loop: do ic = PARAM_FIRST_SCALAR, n_chem
	2232
	2233	CALL horizontal_diffusion_s( chem_tendf(ims,kms,jms,ic), &
	2234	mu, config_flags, &
	2235	chem(ims,kms,jms,ic), &
	2236	msftx, msfty, msfux, msfuy, &
	2237	msfvx, msfvy, xkhh, rdx, rdy, &
	2238	fnm, fnp, cf1, cf2, cf3, &
	2239	zx, zy, rdz, rdzw, dnw, dn, &
	2240	.false., &
	2241	ids, ide, jds, jde, kds, kde, &
	2242	ims, ime, jms, jme, kms, kme, &
	2243	its, ite, jts, jte, kts, kte )
	2244
	2245	ENDDO chem_loop
	2246
	2247	ENDIF
	2248	IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
	2249
	2250	scalar_loop: do is = PARAM_FIRST_SCALAR, n_scalar
	2251
	2252	CALL horizontal_diffusion_s( scalar_tendf(ims,kms,jms,is), &
	2253	mu, config_flags, &
	2254	scalar(ims,kms,jms,is), &
	2255	msftx, msfty, msfux, msfuy, &
	2256	msfvx, msfvy, xkhh, rdx, rdy, &
	2257	fnm, fnp, cf1, cf2, cf3, &
	2258	zx, zy, rdz, rdzw, dnw, dn, &
	2259	.false., &
	2260	ids, ide, jds, jde, kds, kde, &
	2261	ims, ime, jms, jme, kms, kme, &
	2262	its, ite, jts, jte, kts, kte )
	2263
	2264	ENDDO scalar_loop
	2265
	2266	ENDIF
	2267
	2268	END SUBROUTINE horizontal_diffusion_2
	2269
	2270	!=======================================================================
	2271	!=======================================================================
	2272
	2273	SUBROUTINE horizontal_diffusion_u_2( tendency, mu, config_flags, &
	2274	defor11, defor12, div, tke, &
	2275	msfux, msfuy, &
	2276	xkmh, rdx, rdy, fnm, fnp, &
	2277	zx, zy, rdzw, &
	2278	ids, ide, jds, jde, kds, kde, &
	2279	ims, ime, jms, jme, kms, kme, &
	2280	its, ite, jts, jte, kts, kte )
	2281
	2282	!-----------------------------------------------------------------------
	2283	! Begin declarations.
	2284
	2285	IMPLICIT NONE
	2286
	2287	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	2288
	2289	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	2290	ims, ime, jms, jme, kms, kme, &
	2291	its, ite, jts, jte, kts, kte
	2292
	2293	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	2294	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	2295
	2296	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux, &
	2297	msfuy, &
	2298	mu
	2299
	2300	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
	2301
	2302	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdzw
	2303
	2304
	2305	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
	2306	defor12, &
	2307	div, &
	2308	tke, &
	2309	xkmh, &
	2310	zx, &
	2311	zy
	2312
	2313	REAL , INTENT(IN ) :: rdx, &
	2314	rdy
	2315	! Local data
	2316
	2317	INTEGER :: i, j, k, ktf
	2318
	2319	INTEGER :: i_start, i_end, j_start, j_end
	2320	INTEGER :: is_ext,ie_ext,js_ext,je_ext
	2321
	2322	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
	2323	titau2avg, &
	2324	titau1, &
	2325	titau2, &
	2326	xkxavg, &
	2327	rravg
	2328	! new
	2329	! zxavg, &
	2330	! zyavg
	2331	REAL :: mrdx, mrdy, rcoup
	2332
	2333	REAL :: tmpzy, tmpzeta_z
	2334
	2335	REAL :: term1, term2, term3
	2336
	2337	! End declarations.
	2338	!-----------------------------------------------------------------------
	2339
	2340	ktf=MIN(kte,kde-1)
	2341
	2342	!-----------------------------------------------------------------------
	2343	! u : p (.), u(\|), w(-)
	2344	!
	2345	! p u p u u u
	2346	!
	2347	! p \| . \| . \| . \| k+1 \| . \| . \| . \| k+1
	2348	!
	2349	! w - 13 - - k+1 13 k+1
	2350	!
	2351	! p \| 11 O 11 \| . \| k \| 12 O 12 \| . \| k
	2352	!
	2353	! w - 13 - - k 13 k
	2354	!
	2355	! p \| . \| . \| . \| k-1 \| . \| . \| . \| k-1
	2356	!
	2357	! i-1 i i i+1 j-1 j j j+1 j+1
	2358	!
	2359
	2360	i_start = its
	2361	i_end = ite
	2362	j_start = jts
	2363	j_end = MIN(jte,jde-1)
	2364
	2365	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	2366	config_flags%nested) i_start = MAX(ids+1,its)
	2367	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	2368	config_flags%nested) i_end = MIN(ide-1,ite)
	2369	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	2370	config_flags%nested) j_start = MAX(jds+1,jts)
	2371	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	2372	config_flags%nested) j_end = MIN(jde-2,jte)
	2373	IF ( config_flags%periodic_x ) i_start = its
	2374	IF ( config_flags%periodic_x ) i_end = ite
	2375
	2376	! titau1 = titau11
	2377	is_ext=1
	2378	ie_ext=0
	2379	js_ext=0
	2380	je_ext=0
	2381	CALL cal_titau_11_22_33( config_flags, titau1, &
	2382	mu, tke, xkmh, defor11, &
	2383	is_ext, ie_ext, js_ext, je_ext, &
	2384	ids, ide, jds, jde, kds, kde, &
	2385	ims, ime, jms, jme, kms, kme, &
	2386	its, ite, jts, jte, kts, kte )
	2387
	2388	! titau2 = titau12
	2389	is_ext=0
	2390	ie_ext=0
	2391	js_ext=0
	2392	je_ext=1
	2393	CALL cal_titau_12_21( config_flags, titau2, &
	2394	mu, xkmh, defor12, &
	2395	is_ext, ie_ext, js_ext, je_ext, &
	2396	ids, ide, jds, jde, kds, kde, &
	2397	ims, ime, jms, jme, kms, kme, &
	2398	its, ite, jts, jte, kts, kte )
	2399
	2400	! titau1avg = titau11avg
	2401	! titau2avg = titau12avg
	2402
	2403	DO j = j_start, j_end
	2404	DO k = kts+1,ktf
	2405	DO i = i_start, i_end
	2406	titau1avg(i,k,j)=0.5(fnm(k)(titau1(i-1,k ,j)+titau1(i,k ,j))+ &
	2407	fnp(k)*(titau1(i-1,k-1,j)+titau1(i,k-1,j)))
	2408	titau2avg(i,k,j)=0.5(fnm(k)(titau2(i,k ,j+1)+titau2(i,k ,j))+ &
	2409	fnp(k)*(titau2(i,k-1,j+1)+titau2(i,k-1,j)))
	2410	tmpzy = 0.25*( zy(i-1,k,j )+zy(i,k,j )+ &
	2411	zy(i-1,k,j+1)+zy(i,k,j+1) )
	2412	! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
	2413	! titau1avg(i,k,j)=titau1avg(i,k,j)zx(i,k,j)tmpzeta_z
	2414	! titau2avg(i,k,j)=titau2avg(i,k,j)tmpzy tmpzeta_z
	2415
	2416	titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)
	2417	titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
	2418
	2419	ENDDO
	2420	ENDDO
	2421	ENDDO
	2422	!
	2423	DO j = j_start, j_end
	2424	DO i = i_start, i_end
	2425	titau1avg(i,kts,j)=0.
	2426	titau1avg(i,ktf+1,j)=0.
	2427	titau2avg(i,kts,j)=0.
	2428	titau2avg(i,ktf+1,j)=0.
	2429	ENDDO
	2430	ENDDO
	2431	!
	2432	DO j = j_start, j_end
	2433	DO k = kts,ktf
	2434	DO i = i_start, i_end
	2435
	2436	mrdx=msfux(i,j)*rdx
	2437	mrdy=msfuy(i,j)*rdy
	2438	tendency(i,k,j)=tendency(i,k,j)- &
	2439	(mrdx*(titau1(i,k,j )-titau1(i-1,k,j))+ &
	2440	mrdy*(titau2(i,k,j+1)-titau2(i,k,j ))- &
	2441	msfuy(i,j)rdzw(i,k,j)((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
	2442	(titau2avg(i,k+1,j)-titau2avg(i,k,j)) &
	2443	) )
	2444	ENDDO
	2445	ENDDO
	2446	ENDDO
	2447
	2448	END SUBROUTINE horizontal_diffusion_u_2
	2449
	2450	!=======================================================================
	2451	!=======================================================================
	2452
	2453	SUBROUTINE horizontal_diffusion_v_2( tendency, mu, config_flags, &
	2454	defor12, defor22, div, tke, &
	2455	msfvx, msfvy, &
	2456	xkmh, rdx, rdy, fnm, fnp, &
	2457	zx, zy, rdzw, &
	2458	ids, ide, jds, jde, kds, kde, &
	2459	ims, ime, jms, jme, kms, kme, &
	2460	its, ite, jts, jte, kts, kte )
	2461
	2462	!-----------------------------------------------------------------------
	2463	! Begin declarations.
	2464
	2465	IMPLICIT NONE
	2466
	2467	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	2468
	2469	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	2470	ims, ime, jms, jme, kms, kme, &
	2471	its, ite, jts, jte, kts, kte
	2472
	2473	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	2474	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	2475
	2476	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx, &
	2477	msfvy, &
	2478	mu
	2479
	2480	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
	2481
	2482	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor12, &
	2483	defor22, &
	2484	div, &
	2485	tke, &
	2486	xkmh, &
	2487	zx, &
	2488	zy, &
	2489	rdzw
	2490
	2491	REAL , INTENT(IN ) :: rdx, &
	2492	rdy
	2493
	2494	! Local data
	2495
	2496	INTEGER :: i, j, k, ktf
	2497
	2498	INTEGER :: i_start, i_end, j_start, j_end
	2499	INTEGER :: is_ext,ie_ext,js_ext,je_ext
	2500
	2501	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
	2502	titau2avg, &
	2503	titau1, &
	2504	titau2, &
	2505	xkxavg, &
	2506	rravg
	2507	! new
	2508	! zxavg, &
	2509	! zyavg
	2510
	2511	REAL :: mrdx, mrdy, rcoup
	2512
	2513	REAL :: tmpzx, tmpzeta_z
	2514
	2515	! End declarations.
	2516	!-----------------------------------------------------------------------
	2517
	2518	ktf=MIN(kte,kde-1)
	2519
	2520	!-----------------------------------------------------------------------
	2521	! v : p (.), v(+), w(-)
	2522	!
	2523	! p v p v v v
	2524	!
	2525	! p + . + . + . + k+1 + . + . + . + k+1
	2526	!
	2527	! w - 23 - - k+1 23 k+1
	2528	!
	2529	! p + 22 O 22 + . + k + 21 O 21 + . + k
	2530	!
	2531	! w - 23 - - k 23 k
	2532	!
	2533	! p + . + . + . + k-1 + . + . + . + k-1
	2534	!
	2535	! j-1 j j j+1 i-1 i i i+1 i+1
	2536	!
	2537
	2538	i_start = its
	2539	i_end = MIN(ite,ide-1)
	2540	j_start = jts
	2541	j_end = jte
	2542
	2543	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	2544	config_flags%nested) i_start = MAX(ids+1,its)
	2545	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	2546	config_flags%nested) i_end = MIN(ide-2,ite)
	2547	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	2548	config_flags%nested) j_start = MAX(jds+1,jts)
	2549	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	2550	config_flags%nested) j_end = MIN(jde-1,jte)
	2551	IF ( config_flags%periodic_x ) i_start = its
	2552	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
	2553
	2554	! titau1 = titau21
	2555	is_ext=0
	2556	ie_ext=1
	2557	js_ext=0
	2558	je_ext=0
	2559	CALL cal_titau_12_21( config_flags, titau1, &
	2560	mu, xkmh, defor12, &
	2561	is_ext,ie_ext,js_ext,je_ext, &
	2562	ids, ide, jds, jde, kds, kde, &
	2563	ims, ime, jms, jme, kms, kme, &
	2564	its, ite, jts, jte, kts, kte )
	2565
	2566	! titau2 = titau22
	2567	is_ext=0
	2568	ie_ext=0
	2569	js_ext=1
	2570	je_ext=0
	2571	CALL cal_titau_11_22_33( config_flags, titau2, &
	2572	mu, tke, xkmh, defor22, &
	2573	is_ext, ie_ext, js_ext, je_ext, &
	2574	ids, ide, jds, jde, kds, kde, &
	2575	ims, ime, jms, jme, kms, kme, &
	2576	its, ite, jts, jte, kts, kte )
	2577
	2578	DO j = j_start, j_end
	2579	DO k = kts+1,ktf
	2580	DO i = i_start, i_end
	2581	titau1avg(i,k,j)=0.5(fnm(k)(titau1(i+1,k ,j)+titau1(i,k ,j))+ &
	2582	fnp(k)*(titau1(i+1,k-1,j)+titau1(i,k-1,j)))
	2583	titau2avg(i,k,j)=0.5(fnm(k)(titau2(i,k ,j-1)+titau2(i,k ,j))+ &
	2584	fnp(k)*(titau2(i,k-1,j-1)+titau2(i,k-1,j)))
	2585
	2586	tmpzx = 0.25*( zx(i,k,j )+zx(i+1,k,j )+ &
	2587	zx(i,k,j-1)+zx(i+1,k,j-1) )
	2588
	2589
	2590	titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
	2591	titau2avg(i,k,j)=titau2avg(i,k,j)*zy(i,k,j)
	2592
	2593
	2594	ENDDO
	2595	ENDDO
	2596	ENDDO
	2597
	2598	DO j = j_start, j_end
	2599	DO i = i_start, i_end
	2600	titau1avg(i,kts,j)=0.
	2601	titau1avg(i,ktf+1,j)=0.
	2602	titau2avg(i,kts,j)=0.
	2603	titau2avg(i,ktf+1,j)=0.
	2604	ENDDO
	2605	ENDDO
	2606	!
	2607	DO j = j_start, j_end
	2608	DO k = kts,ktf
	2609	DO i = i_start, i_end
	2610
	2611	mrdx=msfvx(i,j)*rdx
	2612	mrdy=msfvy(i,j)*rdy
	2613	tendency(i,k,j)=tendency(i,k,j)- &
	2614	(mrdy*(titau2(i ,k,j)-titau2(i,k,j-1))+ &
	2615	mrdx*(titau1(i+1,k,j)-titau1(i,k,j ))- &
	2616	msfvy(i,j)rdzw(i,k,j)((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
	2617	(titau2avg(i,k+1,j)-titau2avg(i,k,j)) &
	2618	) &
	2619	)
	2620
	2621	ENDDO
	2622	ENDDO
	2623	ENDDO
	2624
	2625	END SUBROUTINE horizontal_diffusion_v_2
	2626
	2627	!=======================================================================
	2628	!=======================================================================
	2629
	2630	SUBROUTINE horizontal_diffusion_w_2( tendency, mu, config_flags, &
	2631	defor13, defor23, div, tke, &
	2632	msftx, msfty, &
	2633	xkmh, rdx, rdy, fnm, fnp, &
	2634	zx, zy, rdz, &
	2635	ids, ide, jds, jde, kds, kde, &
	2636	ims, ime, jms, jme, kms, kme, &
	2637	its, ite, jts, jte, kts, kte )
	2638
	2639	!-----------------------------------------------------------------------
	2640	! Begin declarations.
	2641
	2642	IMPLICIT NONE
	2643
	2644	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	2645
	2646	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	2647	ims, ime, jms, jme, kms, kme, &
	2648	its, ite, jts, jte, kts, kte
	2649
	2650	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	2651	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	2652
	2653	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx, &
	2654	msfty, &
	2655	mu
	2656
	2657	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
	2658
	2659	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
	2660	defor23, &
	2661	div, &
	2662	tke, &
	2663	xkmh, &
	2664	zx, &
	2665	zy, &
	2666	rdz
	2667
	2668	REAL , INTENT(IN ) :: rdx, &
	2669	rdy
	2670
	2671	! Local data
	2672
	2673	INTEGER :: i, j, k, ktf
	2674
	2675	INTEGER :: i_start, i_end, j_start, j_end
	2676	INTEGER :: is_ext,ie_ext,js_ext,je_ext
	2677
	2678	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
	2679	titau2avg, &
	2680	titau1, &
	2681	titau2, &
	2682	xkxavg, &
	2683	rravg
	2684	! new
	2685	! zxavg, &
	2686	! zyavg
	2687
	2688	REAL :: mrdx, mrdy, rcoup
	2689
	2690	REAL :: tmpzx, tmpzy, tmpzeta_z
	2691
	2692	! End declarations.
	2693	!-----------------------------------------------------------------------
	2694
	2695	ktf=MIN(kte,kde-1)
	2696
	2697	!-----------------------------------------------------------------------
	2698	! w : p (.), u(\|), v(+), w(-)
	2699	!
	2700	! p u p u p v p v
	2701	!
	2702	! w - - - k+1 w - - - k+1
	2703	!
	2704	! p . \| 33 \| . k p . + 33 + . k
	2705	!
	2706	! w - 31 O 31 - k w - 32 O 32 - k
	2707	!
	2708	! p . \| 33 \| . k-1 p . \| 33 \| . k-1
	2709	!
	2710	! w - - - k-1 w - - - k-1
	2711	!
	2712	! i-1 i i i+1 j-1 j j j+1
	2713	!
	2714	i_start = its
	2715	i_end = MIN(ite,ide-1)
	2716	j_start = jts
	2717	j_end = MIN(jte,jde-1)
	2718
	2719	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	2720	config_flags%nested) i_start = MAX(ids+1,its)
	2721	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	2722	config_flags%nested) i_end = MIN(ide-2,ite)
	2723	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	2724	config_flags%nested) j_start = MAX(jds+1,jts)
	2725	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	2726	config_flags%nested) j_end = MIN(jde-2,jte)
	2727	IF ( config_flags%periodic_x ) i_start = its
	2728	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
	2729
	2730	! titau1 = titau31
	2731	is_ext=0
	2732	ie_ext=1
	2733	js_ext=0
	2734	je_ext=0
	2735	CALL cal_titau_13_31( config_flags, titau1, defor13, &
	2736	mu, xkmh, fnm, fnp, &
	2737	is_ext, ie_ext, js_ext, je_ext, &
	2738	ids, ide, jds, jde, kds, kde, &
	2739	ims, ime, jms, jme, kms, kme, &
	2740	its, ite, jts, jte, kts, kte )
	2741
	2742	! titau2 = titau32
	2743	is_ext=0
	2744	ie_ext=0
	2745	js_ext=0
	2746	je_ext=1
	2747	CALL cal_titau_23_32( config_flags, titau2, defor23, &
	2748	mu, xkmh, fnm, fnp, &
	2749	is_ext, ie_ext, js_ext, je_ext, &
	2750	ids, ide, jds, jde, kds, kde, &
	2751	ims, ime, jms, jme, kms, kme, &
	2752	its, ite, jts, jte, kts, kte )
	2753
	2754	! titau1avg = titau31avg * zx * zeta_z = titau13avg * zx * zeta_z
	2755	! titau2avg = titau32avg * zy * zeta_z = titau23avg * zy * zeta_z
	2756
	2757	DO j = j_start, j_end
	2758	DO k = kts,ktf
	2759	DO i = i_start, i_end
	2760	titau1avg(i,k,j)=0.25*(titau1(i+1,k+1,j)+titau1(i,k+1,j)+ &
	2761	titau1(i+1,k ,j)+titau1(i,k ,j))
	2762	titau2avg(i,k,j)=0.25*(titau2(i,k+1,j+1)+titau2(i,k+1,j)+ &
	2763	titau2(i,k ,j+1)+titau2(i,k ,j))
	2764	! new
	2765	tmpzx =0.25*( zx(i,k ,j)+zx(i+1,k ,j)+ &
	2766	zx(i,k+1,j)+zx(i+1,k+1,j) )
	2767	tmpzy =0.25*( zy(i,k ,j)+zy(i,k ,j+1)+ &
	2768	zy(i,k+1,j)+zy(i,k+1,j+1) )
	2769
	2770	titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
	2771	titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
	2772	! titau1avg(i,k,j)=titau1avg(i,k,j)tmpzxzeta_z(i,j)
	2773	! titau2avg(i,k,j)=titau2avg(i,k,j)tmpzyzeta_z(i,j)
	2774	ENDDO
	2775	ENDDO
	2776	ENDDO
	2777
	2778	DO j = j_start, j_end
	2779	DO i = i_start, i_end
	2780	titau1avg(i,ktf+1,j)=0.
	2781	titau2avg(i,ktf+1,j)=0.
	2782	ENDDO
	2783	ENDDO
	2784
	2785	DO j = j_start, j_end
	2786	DO k = kts+1,ktf
	2787	DO i = i_start, i_end
	2788
	2789	mrdx=msftx(i,j)*rdx
	2790	mrdy=msfty(i,j)*rdy
	2791
	2792	tendency(i,k,j)=tendency(i,k,j)- &
	2793	(mrdx*(titau1(i+1,k,j)-titau1(i,k,j))+ &
	2794	mrdy*(titau2(i,k,j+1)-titau2(i,k,j))- &
	2795	msfty(i,j)rdz(i,k,j)(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
	2796	titau2avg(i,k,j)-titau2avg(i,k-1,j) &
	2797	) &
	2798	)
	2799	! msft(i,j)/dn(k)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
	2800	! titau2avg(i,k,j)-titau2avg(i,k-1,j) &
	2801	! ) &
	2802	! )
	2803	ENDDO
	2804	ENDDO
	2805	ENDDO
	2806
	2807	END SUBROUTINE horizontal_diffusion_w_2
	2808
	2809	!=======================================================================
	2810	!=======================================================================
	2811
	2812	SUBROUTINE horizontal_diffusion_s (tendency, mu, config_flags, var, &
	2813	msftx, msfty, msfux, msfuy, &
	2814	msfvx, msfvy, xkhh, rdx, rdy, &
	2815	fnm, fnp, cf1, cf2, cf3, &
	2816	zx, zy, rdz, rdzw, dnw, dn, &
	2817	doing_tke, &
	2818	ids, ide, jds, jde, kds, kde, &
	2819	ims, ime, jms, jme, kms, kme, &
	2820	its, ite, jts, jte, kts, kte )
	2821
	2822	!-----------------------------------------------------------------------
	2823	! Begin declarations.
	2824
	2825	IMPLICIT NONE
	2826
	2827	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	2828
	2829	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	2830	ims, ime, jms, jme, kms, kme, &
	2831	its, ite, jts, jte, kts, kte
	2832
	2833	LOGICAL, INTENT(IN ) :: doing_tke
	2834
	2835	REAL , INTENT(IN ) :: cf1, cf2, cf3
	2836
	2837	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	2838	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	2839	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
	2840	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
	2841
	2842	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux
	2843	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfuy
	2844	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx
	2845	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvy
	2846	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx
	2847	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfty
	2848
	2849	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: mu
	2850
	2851	! REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
	2852	! INTENT(IN ) :: xkhh
	2853
	2854	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
	2855
	2856	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: &
	2857	xkhh, &
	2858	rdz, &
	2859	rdzw
	2860
	2861	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: var, &
	2862	zx, &
	2863	zy
	2864
	2865	REAL , INTENT(IN ) :: rdx, &
	2866	rdy
	2867
	2868	! Local data
	2869
	2870	INTEGER :: i, j, k, ktf
	2871
	2872	INTEGER :: i_start, i_end, j_start, j_end
	2873
	2874	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: H1avg, &
	2875	H2avg, &
	2876	H1, &
	2877	H2, &
	2878	xkxavg
	2879	! new
	2880	! zxavg, &
	2881	! zyavg
	2882
	2883	REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
	2884
	2885	REAL :: mrdx, mrdy, rcoup
	2886	REAL :: tmpzx, tmpzy, tmpzeta_z, rdzu, rdzv
	2887	INTEGER :: ktes1,ktes2
	2888
	2889	! End declarations.
	2890	!-----------------------------------------------------------------------
	2891
	2892	ktf=MIN(kte,kde-1)
	2893
	2894	!-----------------------------------------------------------------------
	2895	! scalars: t (.), u(\|), v(+), w(-)
	2896	!
	2897	! t u t u t v t v
	2898	!
	2899	! w - 3 - k+1 w - 3 - k+1
	2900	!
	2901	! t . 1 O 1 . k t . 2 O 2 . k
	2902	!
	2903	! w - 3 - k w - 3 - k
	2904	!
	2905	! t . \| . \| . k-1 t . + . + . k-1
	2906	!
	2907	! w - - - k-1 w - - - k-1
	2908	!
	2909	! t i-1 i i i+1 j-1 j j j+1
	2910	!
	2911
	2912	ktes1=kte-1
	2913	ktes2=kte-2
	2914
	2915	i_start = its
	2916	i_end = MIN(ite,ide-1)
	2917	j_start = jts
	2918	j_end = MIN(jte,jde-1)
	2919
	2920	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	2921	config_flags%nested) i_start = MAX(ids+1,its)
	2922	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	2923	config_flags%nested) i_end = MIN(ide-2,ite)
	2924	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	2925	config_flags%nested) j_start = MAX(jds+1,jts)
	2926	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	2927	config_flags%nested) j_end = MIN(jde-2,jte)
	2928	IF ( config_flags%periodic_x ) i_start = its
	2929	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
	2930
	2931	! diffusion of the TKE needs mutiple 2
	2932
	2933	IF ( doing_tke ) THEN
	2934	DO j = j_start, j_end
	2935	DO k = kts,ktf
	2936	DO i = i_start, i_end
	2937	tmptendf(i,k,j)=tendency(i,k,j)
	2938	ENDDO
	2939	ENDDO
	2940	ENDDO
	2941	ENDIF
	2942
	2943	! H1 = partial var over partial x
	2944
	2945	DO j = j_start, j_end
	2946	DO k = kts, ktf
	2947	DO i = i_start, i_end + 1
	2948	! new
	2949	! zxavg(i,k,j) =0.5*( zx(i-1,k,j)+ zx(i,k,j))
	2950	xkxavg(i,k,j)=0.5*(xkhh(i-1,k,j)+xkhh(i,k,j))
	2951	ENDDO
	2952	ENDDO
	2953	ENDDO
	2954
	2955	DO j = j_start, j_end
	2956	DO k = kts+1, ktf
	2957	DO i = i_start, i_end + 1
	2958	H1avg(i,k,j)=0.5(fnm(k)(var(i-1,k ,j)+var(i,k ,j))+ &
	2959	fnp(k)*(var(i-1,k-1,j)+var(i,k-1,j)))
	2960	ENDDO
	2961	ENDDO
	2962	ENDDO
	2963
	2964	DO j = j_start, j_end
	2965	DO i = i_start, i_end + 1
	2966	H1avg(i,kts ,j)=0.5(cf1var(i ,1,j)+cf2*var(i ,2,j)+ &
	2967	cf3var(i ,3,j)+cf1var(i-1,1,j)+ &
	2968	cf2var(i-1,2,j)+cf3var(i-1,3,j))
	2969	H1avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
	2970	var(i,ktes2,j))0.5dnw(ktes1)/dn(ktes1)+ &
	2971	var(i-1,ktes1,j)+(var(i-1,ktes1,j)- &
	2972	var(i-1,ktes2,j))0.5dnw(ktes1)/dn(ktes1))
	2973	ENDDO
	2974	ENDDO
	2975
	2976	DO j = j_start, j_end
	2977	DO k = kts, ktf
	2978	DO i = i_start, i_end + 1
	2979	! new
	2980	tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j))
	2981	rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
	2982	H1(i,k,j)=-msfuy(i,j)xkxavg(i,k,j)( &
	2983	rdx(var(i,k,j)-var(i-1,k,j)) - tmpzx &
	2984	(H1avg(i,k+1,j)-H1avg(i,k,j))*rdzu )
	2985
	2986	! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
	2987	! H1(i,k,j)=-msfuy(i,j)xkxavg(i,k,j)( &
	2988	! rdx(var(i,k,j)-var(i-1,k,j)) - tmpzxtmpzeta_z* &
	2989	! (H1avg(i,k+1,j)-H1avg(i,k,j))/dnw(k))
	2990	ENDDO
	2991	ENDDO
	2992	ENDDO
	2993
	2994	! H2 = partial var over partial y
	2995
	2996	DO j = j_start, j_end + 1
	2997	DO k = kts, ktf
	2998	DO i = i_start, i_end
	2999	! new
	3000	! zyavg(i,k,j) =0.5*( zy(i,k,j-1)+ zy(i,k,j))
	3001	xkxavg(i,k,j)=0.5*(xkhh(i,k,j-1)+xkhh(i,k,j))
	3002	ENDDO
	3003	ENDDO
	3004	ENDDO
	3005
	3006	DO j = j_start, j_end + 1
	3007	DO k = kts+1, ktf
	3008	DO i = i_start, i_end
	3009	! new
	3010	H2avg(i,k,j)=0.5(fnm(k)(var(i,k ,j-1)+var(i,k ,j))+ &
	3011	fnp(k)*(var(i,k-1,j-1)+var(i,k-1,j)))
	3012	ENDDO
	3013	ENDDO
	3014	ENDDO
	3015
	3016	DO j = j_start, j_end + 1
	3017	DO i = i_start, i_end
	3018	H2avg(i,kts ,j)=0.5(cf1var(i,1,j )+cf2*var(i ,2,j)+ &
	3019	cf3var(i,3,j )+cf1var(i,1,j-1)+ &
	3020	cf2var(i,2,j-1)+cf3var(i,3,j-1))
	3021	H2avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
	3022	var(i,ktes2,j))0.5dnw(ktes1)/dn(ktes1)+ &
	3023	var(i,ktes1,j-1)+(var(i,ktes1,j-1)- &
	3024	var(i,ktes2,j-1))0.5dnw(ktes1)/dn(ktes1))
	3025	ENDDO
	3026	ENDDO
	3027
	3028	DO j = j_start, j_end + 1
	3029	DO k = kts, ktf
	3030	DO i = i_start, i_end
	3031	! new
	3032	tmpzy = 0.5*( zy(i,k,j)+ zy(i,k+1,j))
	3033	rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
	3034	H2(i,k,j)=-msfvy(i,j)xkxavg(i,k,j)( &
	3035	rdy(var(i,k,j)-var(i,k,j-1)) - tmpzy &
	3036	(H2avg(i ,k+1,j)-H2avg(i,k,j))*rdzv)
	3037
	3038	! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i,j-1))
	3039	! H2(i,k,j)=-msfvy(i,j)xkxavg(i,k,j)( &
	3040	! rdy(var(i,k,j)-var(i,k,j-1)) - tmpzytmpzeta_z* &
	3041	! (H2avg(i ,k+1,j)-H2avg(i,k,j))/dnw(k))
	3042	ENDDO
	3043	ENDDO
	3044	ENDDO
	3045
	3046	DO j = j_start, j_end
	3047	DO k = kts+1, ktf
	3048	DO i = i_start, i_end
	3049	H1avg(i,k,j)=0.5(fnm(k)(H1(i+1,k ,j)+H1(i,k ,j))+ &
	3050	fnp(k)*(H1(i+1,k-1,j)+H1(i,k-1,j)))
	3051	H2avg(i,k,j)=0.5(fnm(k)(H2(i,k ,j+1)+H2(i,k ,j))+ &
	3052	fnp(k)*(H2(i,k-1,j+1)+H2(i,k-1,j)))
	3053	! new
	3054	! zxavg(i,k,j)=fnm(k)zx(i,k,j)+fnp(k)zx(i,k-1,j)
	3055	! zyavg(i,k,j)=fnm(k)zy(i,k,j)+fnp(k)zy(i,k-1,j)
	3056
	3057	! H1avg(i,k,j)=zxH1avgzeta_z
	3058	! H2avg(i,k,j)=zyH2avgzeta_z
	3059
	3060	tmpzx = 0.5*( zx(i,k,j)+ zx(i+1,k,j ))
	3061	tmpzy = 0.5*( zy(i,k,j)+ zy(i ,k,j+1))
	3062
	3063	H1avg(i,k,j)=H1avg(i,k,j)*tmpzx
	3064	H2avg(i,k,j)=H2avg(i,k,j)*tmpzy
	3065
	3066	! H1avg(i,k,j)=H1avg(i,k,j)tmpzxzeta_z(i,j)
	3067	! H2avg(i,k,j)=H2avg(i,k,j)tmpzyzeta_z(i,j)
	3068	ENDDO
	3069	ENDDO
	3070	ENDDO
	3071
	3072	DO j = j_start, j_end
	3073	DO i = i_start, i_end
	3074	H1avg(i,kts ,j)=0.
	3075	H1avg(i,ktf+1,j)=0.
	3076	H2avg(i,kts ,j)=0.
	3077	H2avg(i,ktf+1,j)=0.
	3078	ENDDO
	3079	ENDDO
	3080
	3081	DO j = j_start, j_end
	3082	DO k = kts,ktf
	3083	DO i = i_start, i_end
	3084
	3085	mrdx=msftx(i,j)*rdx
	3086	mrdy=msfty(i,j)*rdy
	3087
	3088	tendency(i,k,j)=tendency(i,k,j)- &
	3089	(mrdx0.5((mu(i+1,j)+mu(i,j))*H1(i+1,k,j)- &
	3090	(mu(i-1,j)+mu(i,j))*H1(i ,k,j))+ &
	3091	mrdy0.5((mu(i,j+1)+mu(i,j))*H2(i,k,j+1)- &
	3092	(mu(i,j-1)+mu(i,j))*H2(i,k,j ))- &
	3093	msfty(i,j)mu(i,j)(H1avg(i,k+1,j)-H1avg(i,k,j)+ &
	3094	H2avg(i,k+1,j)-H2avg(i,k,j) &
	3095	)*rdzw(i,k,j) &
	3096	)
	3097
	3098	ENDDO
	3099	ENDDO
	3100	ENDDO
	3101
	3102	IF ( doing_tke ) THEN
	3103	DO j = j_start, j_end
	3104	DO k = kts,ktf
	3105	DO i = i_start, i_end
	3106	tendency(i,k,j)=tmptendf(i,k,j)+2.* &
	3107	(tendency(i,k,j)-tmptendf(i,k,j))
	3108	ENDDO
	3109	ENDDO
	3110	ENDDO
	3111	ENDIF
	3112
	3113	END SUBROUTINE horizontal_diffusion_s
	3114
	3115	!=======================================================================
	3116	!=======================================================================
	3117
	3118	SUBROUTINE vertical_diffusion_2 ( ru_tendf, rv_tendf, rw_tendf, rt_tendf, &
	3119	tke_tendf, moist_tendf, n_moist, &
	3120	chem_tendf, n_chem, &
	3121	scalar_tendf, n_scalar, &
	3122	u_2, v_2, &
	3123	thp,u_base,v_base,t_base,qv_base,mu,tke, &
	3124	config_flags,defor13,defor23,defor33,div, &
	3125	moist, chem, scalar, xkmv, xkhv,km_opt, &
	3126	fnm, fnp, dn, dnw, rdz, rdzw, &
	3127	hfx, qfx, ust, rho, &
	3128	ids, ide, jds, jde, kds, kde, &
	3129	ims, ime, jms, jme, kms, kme, &
	3130	its, ite, jts, jte, kts, kte )
	3131
	3132	!-----------------------------------------------------------------------
	3133	! Begin declarations.
	3134
	3135	IMPLICIT NONE
	3136
	3137	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	3138
	3139	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	3140	ims, ime, jms, jme, kms, kme, &
	3141	its, ite, jts, jte, kts, kte
	3142
	3143	INTEGER , INTENT(IN ) :: n_moist, n_chem, n_scalar, km_opt
	3144
	3145	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	3146	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	3147	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
	3148	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
	3149	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
	3150
	3151	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: qv_base
	3152	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base
	3153	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: v_base
	3154	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: t_base
	3155
	3156	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::ru_tendf,&
	3157	rv_tendf,&
	3158	rw_tendf,&
	3159	tke_tendf,&
	3160	rt_tendf
	3161
	3162	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
	3163	INTENT(INOUT) :: moist_tendf
	3164
	3165	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
	3166	INTENT(INOUT) :: chem_tendf
	3167
	3168	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
	3169	INTENT(INOUT) :: scalar_tendf
	3170
	3171	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
	3172	INTENT(INOUT) :: moist
	3173
	3174	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
	3175	INTENT(INOUT) :: chem
	3176
	3177	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
	3178	INTENT(IN ) :: scalar
	3179
	3180	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
	3181	defor23, &
	3182	defor33, &
	3183	div, &
	3184	xkmv, &
	3185	xkhv, &
	3186	tke, &
	3187	rdz, &
	3188	u_2, &
	3189	v_2, &
	3190	rdzw
	3191
	3192	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rho
	3193	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: hfx, &
	3194	qfx
	3195	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: ust
	3196	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: thp
	3197
	3198	! LOCAL VAR
	3199
	3200	REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: var_mix
	3201
	3202	INTEGER :: im, i,j,k
	3203	INTEGER :: i_start, i_end, j_start, j_end
	3204
	3205	! REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: xkhv
	3206
	3207	!***************************************************************************
	3208	!***************************************************************************
	3209	!MODIFICA VARIABILI PER I FLUSSI
	3210	!
	3211	REAL :: V0_u,V0_v,tao_xz,tao_yz,ustar,cd0
	3212	REAL :: xsfc,psi1,vk2,zrough,lnz
	3213	REAL :: heat_flux, moist_flux, heat_flux0
	3214	REAL :: cpm
	3215	!
	3216	!FINE MODIFICA VARIABILI PER I FLUSSI
	3217	!***************************************************************************
	3218	!
	3219
	3220	! End declarations.
	3221	!-----------------------------------------------------------------------
	3222
	3223	i_start = its
	3224	i_end = MIN(ite,ide-1)
	3225	j_start = jts
	3226	j_end = MIN(jte,jde-1)
	3227	!
	3228	!-----------------------------------------------------------------------
	3229
	3230	CALL vertical_diffusion_u_2( ru_tendf, config_flags, mu, &
	3231	defor13, xkmv, &
	3232	dnw, rdzw, fnm, fnp, &
	3233	ids, ide, jds, jde, kds, kde, &
	3234	ims, ime, jms, jme, kms, kme, &
	3235	its, ite, jts, jte, kts, kte )
	3236
	3237
	3238	CALL vertical_diffusion_v_2( rv_tendf, config_flags, mu, &
	3239	defor23, xkmv, &
	3240	dnw, rdzw, fnm, fnp, &
	3241	ids, ide, jds, jde, kds, kde, &
	3242	ims, ime, jms, jme, kms, kme, &
	3243	its, ite, jts, jte, kts, kte )
	3244
	3245	CALL vertical_diffusion_w_2( rw_tendf, config_flags, mu, &
	3246	defor33, tke(ims,kms,jms), &
	3247	div, xkmv, &
	3248	dn, rdz, &
	3249	ids, ide, jds, jde, kds, kde, &
	3250	ims, ime, jms, jme, kms, kme, &
	3251	its, ite, jts, jte, kts, kte )
	3252
	3253	!*****************************************
	3254	!*****************************************
	3255	! MODIFICA al flusso di momento alla parete
	3256	!
	3257	vflux: SELECT CASE( config_flags%isfflx )
	3258	CASE (0) ! Assume cd a constant, specified in namelist
	3259	cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient
	3260	! set in namelist.input
	3261	!
	3262	!calcolo del modulo della velocita
	3263	DO j = j_start, j_end
	3264	DO i = i_start, ite
	3265	V0_u=0.
	3266	tao_xz=0.
	3267	V0_u= sqrt((u_2(i,kts,j)**2) + &
	3268	(((v_2(i ,kts,j )+ &
	3269	v_2(i ,kts,j+1)+ &
	3270	v_2(i-1,kts,j )+ &
	3271	v_2(i-1,kts,j+1))/4)**2))+epsilon
	3272	tao_xz=cd0V0_uu_2(i,kts,j)
	3273	ru_tendf(i,kts,j)=ru_tendf(i,kts,j) &
	3274	-0.25(mu(i,j)+mu(i-1,j))tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
	3275	ENDDO
	3276	ENDDO
	3277
	3278	!
	3279	DO j = j_start, jte
	3280	DO i = i_start, i_end
	3281	V0_v=0.
	3282	tao_yz=0.
	3283	V0_v= sqrt((v_2(i,kts,j)**2) + &
	3284	(((u_2(i ,kts,j )+ &
	3285	u_2(i ,kts,j-1)+ &
	3286	u_2(i+1,kts,j )+ &
	3287	u_2(i+1,kts,j-1))/4)**2))+epsilon
	3288	tao_yz=cd0V0_vv_2(i,kts,j)
	3289	rv_tendf(i,kts,j)=rv_tendf(i,kts,j) &
	3290	-0.25(mu(i,j)+mu(i,j-1))tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
	3291	ENDDO
	3292	ENDDO
	3293
	3294	CASE (1,2) ! ustar computed from surface routine
	3295	DO j = j_start, j_end
	3296	DO i = i_start, ite
	3297	V0_u=0.
	3298	tao_xz=0.
	3299	V0_u= sqrt((u_2(i,kts,j)**2) + &
	3300	(((v_2(i ,kts,j )+ &
	3301	v_2(i ,kts,j+1)+ &
	3302	v_2(i-1,kts,j )+ &
	3303	v_2(i-1,kts,j+1))/4)**2))+epsilon
	3304	ustar=0.5*(ust(i,j)+ust(i-1,j))
	3305	tao_xz=ustarustaru_2(i,kts,j)/V0_u
	3306	ru_tendf(i,kts,j)=ru_tendf(i,kts,j) &
	3307	-0.25(mu(i,j)+mu(i-1,j))tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
	3308	ENDDO
	3309	ENDDO
	3310
	3311	DO j = j_start, jte
	3312	DO i = i_start, i_end
	3313	V0_v=0.
	3314	tao_yz=0.
	3315	V0_v= sqrt((v_2(i,kts,j)**2) + &
	3316	(((u_2(i ,kts,j )+ &
	3317	u_2(i ,kts,j-1)+ &
	3318	u_2(i+1,kts,j )+ &
	3319	u_2(i+1,kts,j-1))/4)**2))+epsilon
	3320	ustar=0.5*(ust(i,j)+ust(i,j-1))
	3321	tao_yz=ustarustarv_2(i,kts,j)/V0_v
	3322	rv_tendf(i,kts,j)=rv_tendf(i,kts,j) &
	3323	-0.25(mu(i,j)+mu(i,j-1))tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
	3324	ENDDO
	3325	ENDDO
	3326
	3327	CASE DEFAULT
	3328	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
	3329	END SELECT vflux
	3330
	3331	!
	3332	! FINE MODIFICA al flusso di momento alla parete
	3333	!*****************************************
	3334	!*****************************************
	3335
	3336	IF ( config_flags%mix_full_fields ) THEN
	3337
	3338	DO j=jts,min(jte,jde-1)
	3339	DO k=kts,kte-1
	3340	DO i=its,min(ite,ide-1)
	3341	var_mix(i,k,j) = thp(i,k,j)
	3342	ENDDO
	3343	ENDDO
	3344	ENDDO
	3345
	3346	ELSE
	3347
	3348	DO j=jts,min(jte,jde-1)
	3349	DO k=kts,kte-1
	3350	DO i=its,min(ite,ide-1)
	3351	var_mix(i,k,j) = thp(i,k,j) - t_base(k)
	3352	ENDDO
	3353	ENDDO
	3354	ENDDO
	3355
	3356	END IF
	3357
	3358	CALL vertical_diffusion_s( rt_tendf, config_flags, var_mix, mu, xkhv, &
	3359	dn, dnw, rdz, rdzw, fnm, fnp, &
	3360	.false., &
	3361	ids, ide, jds, jde, kds, kde, &
	3362	ims, ime, jms, jme, kms, kme, &
	3363	its, ite, jts, jte, kts, kte )
	3364
	3365
	3366	!*****************************************
	3367	!*****************************************
	3368	!MODIFICA al flusso di calore
	3369	!
	3370	!
	3371	hflux: SELECT CASE( config_flags%isfflx )
	3372	CASE (0,2) ! with fixed surface heat flux given in the namelist
	3373	heat_flux = config_flags%tke_heat_flux ! constant heat flux value
	3374	! set in namelist.input
	3375	DO j = j_start, j_end
	3376	DO i = i_start, i_end
	3377	rt_tendf(i,kts,j)=rt_tendf(i,kts,j) &
	3378	+mu(i,j)heat_fluxrdzw(i,kts,j)
	3379	ENDDO
	3380	ENDDO
	3381
	3382	CASE (1) ! use surface heat flux computed from surface routine
	3383	DO j = j_start, j_end
	3384	DO i = i_start, i_end
	3385
	3386	cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV))
	3387	heat_flux = hfx(i,j)/cpm/rho(i,1,j)
	3388	rt_tendf(i,kts,j)=rt_tendf(i,kts,j) &
	3389	+mu(i,j)heat_fluxrdzw(i,kts,j)
	3390
	3391	ENDDO
	3392	ENDDO
	3393
	3394	CASE DEFAULT
	3395	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
	3396	END SELECT hflux
	3397
	3398	!
	3399	! FINE MODIFICA al flusso di calore
	3400	!*****************************************
	3401	!*****************************************
	3402
	3403	If (km_opt .eq. 2) then
	3404	CALL vertical_diffusion_s( tke_tendf(ims,kms,jms), &
	3405	config_flags, tke(ims,kms,jms), &
	3406	mu, xkhv, &
	3407	dn, dnw, rdz, rdzw, fnm, fnp, &
	3408	.true., &
	3409	ids, ide, jds, jde, kds, kde, &
	3410	ims, ime, jms, jme, kms, kme, &
	3411	its, ite, jts, jte, kts, kte )
	3412	endif
	3413
	3414	IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
	3415
	3416	moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
	3417
	3418	IF ( (.not. config_flags%mix_full_fields) .and. (im == P_QV) ) THEN
	3419
	3420	DO j=jts,min(jte,jde-1)
	3421	DO k=kts,kte-1
	3422	DO i=its,min(ite,ide-1)
	3423	var_mix(i,k,j) = moist(i,k,j,im) - qv_base(k)
	3424	ENDDO
	3425	ENDDO
	3426	ENDDO
	3427
	3428	ELSE
	3429
	3430	DO j=jts,min(jte,jde-1)
	3431	DO k=kts,kte-1
	3432	DO i=its,min(ite,ide-1)
	3433	var_mix(i,k,j) = moist(i,k,j,im)
	3434	ENDDO
	3435	ENDDO
	3436	ENDDO
	3437
	3438	END IF
	3439
	3440
	3441	CALL vertical_diffusion_s( moist_tendf(ims,kms,jms,im), &
	3442	config_flags, var_mix, &
	3443	mu, xkhv, &
	3444	dn, dnw, rdz, rdzw, fnm, fnp, &
	3445	.false., &
	3446	ids, ide, jds, jde, kds, kde, &
	3447	ims, ime, jms, jme, kms, kme, &
	3448	its, ite, jts, jte, kts, kte )
	3449
	3450	!*****************************************
	3451	!*****************************************
	3452	!MODIFICATIONS for water vapor flux
	3453	!
	3454	!
	3455	qflux: SELECT CASE( config_flags%isfflx )
	3456	CASE (0)
	3457	! do nothing
	3458	CASE (1,2) ! with surface moisture flux
	3459	IF ( im == P_QV ) THEN
	3460	DO j = j_start, j_end
	3461	DO i = i_start, i_end
	3462	moist_flux = qfx(i,j)/rho(i,1,j)
	3463	moist_tendf(i,kts,j,im)=moist_tendf(i,kts,j,im) &
	3464	+mu(i,j)moist_fluxrdzw(i,kts,j)
	3465	ENDDO
	3466	ENDDO
	3467	ENDIF
	3468
	3469	CASE DEFAULT
	3470	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
	3471	END SELECT qflux
	3472	!
	3473	! END MODIFICATIONS for water vapor flux
	3474	!*****************************************
	3475	!*****************************************
	3476	ENDDO moist_loop
	3477
	3478	ENDIF
	3479
	3480	IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN
	3481
	3482	chem_loop: do im = PARAM_FIRST_SCALAR, n_chem
	3483
	3484	CALL vertical_diffusion_s( chem_tendf(ims,kms,jms,im), &
	3485	config_flags, chem(ims,kms,jms,im), &
	3486	mu, xkhv, &
	3487	dn, dnw, rdz, rdzw, fnm, fnp, &
	3488	.false., &
	3489	ids, ide, jds, jde, kds, kde, &
	3490	ims, ime, jms, jme, kms, kme, &
	3491	its, ite, jts, jte, kts, kte )
	3492	ENDDO chem_loop
	3493
	3494	ENDIF
	3495
	3496
	3497	IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
	3498
	3499	scalar_loop: do im = PARAM_FIRST_SCALAR, n_scalar
	3500
	3501	CALL vertical_diffusion_s( scalar_tendf(ims,kms,jms,im), &
	3502	config_flags, scalar(ims,kms,jms,im), &
	3503	mu, xkhv, &
	3504	dn, dnw, rdz, rdzw, fnm, fnp, &
	3505	.false., &
	3506	ids, ide, jds, jde, kds, kde, &
	3507	ims, ime, jms, jme, kms, kme, &
	3508	its, ite, jts, jte, kts, kte )
	3509	ENDDO scalar_loop
	3510
	3511	ENDIF
	3512
	3513	END SUBROUTINE vertical_diffusion_2
	3514
	3515	!=======================================================================
	3516	!=======================================================================
	3517
	3518	SUBROUTINE vertical_diffusion_u_2( tendency, config_flags, mu, &
	3519	defor13, xkmv, &
	3520	dnw, rdzw, fnm, fnp, &
	3521	ids, ide, jds, jde, kds, kde, &
	3522	ims, ime, jms, jme, kms, kme, &
	3523	its, ite, jts, jte, kts, kte )
	3524
	3525	!-----------------------------------------------------------------------
	3526	! Begin declarations.
	3527
	3528	IMPLICIT NONE
	3529
	3530	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	3531
	3532	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	3533	ims, ime, jms, jme, kms, kme, &
	3534	its, ite, jts, jte, kts, kte
	3535
	3536	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	3537	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	3538	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
	3539	! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
	3540
	3541	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
	3542
	3543	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
	3544	INTENT(IN ) ::defor13, &
	3545	xkmv, &
	3546	rdzw
	3547	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
	3548
	3549	! LOCAL VARS
	3550
	3551	INTEGER :: i, j, k, ktf
	3552
	3553	INTEGER :: i_start, i_end, j_start, j_end
	3554	INTEGER :: is_ext,ie_ext,js_ext,je_ext
	3555
	3556	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
	3557
	3558	REAL , DIMENSION( its:ite, jts:jte) :: zzavg
	3559
	3560	REAL :: rdzu
	3561
	3562	! End declarations.
	3563	!-----------------------------------------------------------------------
	3564
	3565	ktf=MIN(kte,kde-1)
	3566
	3567	i_start = its
	3568	i_end = ite
	3569	j_start = jts
	3570	j_end = MIN(jte,jde-1)
	3571
	3572	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	3573	config_flags%nested) i_start = MAX(ids+1,its)
	3574	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	3575	config_flags%nested) i_end = MIN(ide-1,ite)
	3576	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	3577	config_flags%nested) j_start = MAX(jds+1,jts)
	3578	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	3579	config_flags%nested) j_end = MIN(jde-2,jte)
	3580	IF ( config_flags%periodic_x ) i_start = its
	3581	IF ( config_flags%periodic_x ) i_end = ite
	3582
	3583	! titau3 = titau13
	3584	is_ext=0
	3585	ie_ext=0
	3586	js_ext=0
	3587	je_ext=0
	3588	CALL cal_titau_13_31( config_flags, titau3, defor13, &
	3589	mu, xkmv, fnm, fnp, &
	3590	is_ext, ie_ext, js_ext, je_ext, &
	3591	ids, ide, jds, jde, kds, kde, &
	3592	ims, ime, jms, jme, kms, kme, &
	3593	its, ite, jts, jte, kts, kte )
	3594	!
	3595	DO j = j_start, j_end
	3596	DO k=kts+1,ktf
	3597	DO i = i_start, i_end
	3598
	3599	rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
	3600	tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j)-titau3(i,k,j))
	3601
	3602	ENDDO
	3603	ENDDO
	3604	ENDDO
	3605
	3606	! ******** MODIF...
	3607	! we will pick up the surface drag (titau3(i,kts,j)) later
	3608	!
	3609	DO j = j_start, j_end
	3610	k=kts
	3611	DO i = i_start, i_end
	3612
	3613	rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
	3614	tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j))
	3615	ENDDO
	3616	ENDDO
	3617	! ******** MODIF...
	3618
	3619	END SUBROUTINE vertical_diffusion_u_2
	3620
	3621	!=======================================================================
	3622	!=======================================================================
	3623
	3624	SUBROUTINE vertical_diffusion_v_2( tendency, config_flags, mu, &
	3625	defor23, xkmv, &
	3626	dnw, rdzw, fnm, fnp, &
	3627	ids, ide, jds, jde, kds, kde, &
	3628	ims, ime, jms, jme, kms, kme, &
	3629	its, ite, jts, jte, kts, kte )
	3630	!-----------------------------------------------------------------------
	3631	! Begin declarations.
	3632
	3633	IMPLICIT NONE
	3634
	3635	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	3636
	3637	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	3638	ims, ime, jms, jme, kms, kme, &
	3639	its, ite, jts, jte, kts, kte
	3640	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	3641	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	3642	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
	3643	! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
	3644
	3645	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
	3646
	3647	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
	3648	INTENT(IN ) ::defor23, &
	3649	xkmv, &
	3650	rdzw
	3651
	3652	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
	3653
	3654	! LOCAL VARS
	3655
	3656	INTEGER :: i, j, k, ktf
	3657
	3658	INTEGER :: i_start, i_end, j_start, j_end
	3659	INTEGER :: is_ext,ie_ext,js_ext,je_ext
	3660
	3661	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
	3662
	3663	REAL , DIMENSION( its:ite, jts:jte) :: zzavg
	3664
	3665	REAL :: rdzv
	3666
	3667	! End declarations.
	3668	!-----------------------------------------------------------------------
	3669
	3670	ktf=MIN(kte,kde-1)
	3671
	3672	i_start = its
	3673	i_end = MIN(ite,ide-1)
	3674	j_start = jts
	3675	j_end = jte
	3676
	3677	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	3678	config_flags%nested) i_start = MAX(ids+1,its)
	3679	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	3680	config_flags%nested) i_end = MIN(ide-2,ite)
	3681	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	3682	config_flags%nested) j_start = MAX(jds+1,jts)
	3683	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	3684	config_flags%nested) j_end = MIN(jde-1,jte)
	3685	IF ( config_flags%periodic_x ) i_start = its
	3686	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
	3687
	3688	! titau3 = titau23
	3689	is_ext=0
	3690	ie_ext=0
	3691	js_ext=0
	3692	je_ext=0
	3693	CALL cal_titau_23_32( config_flags, titau3, defor23, &
	3694	mu, xkmv, fnm, fnp, &
	3695	is_ext, ie_ext, js_ext, je_ext, &
	3696	ids, ide, jds, jde, kds, kde, &
	3697	ims, ime, jms, jme, kms, kme, &
	3698	its, ite, jts, jte, kts, kte )
	3699
	3700	DO j = j_start, j_end
	3701	DO k = kts+1,ktf
	3702	DO i = i_start, i_end
	3703
	3704	rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
	3705	tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j)-titau3(i,k,j))
	3706
	3707	ENDDO
	3708	ENDDO
	3709	ENDDO
	3710
	3711	! ******** MODIF...
	3712	! we will pick up the surface drag (titau3(i,kts,j)) later
	3713	!
	3714	DO j = j_start, j_end
	3715	k=kts
	3716	DO i = i_start, i_end
	3717
	3718	rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
	3719	tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j))
	3720
	3721	ENDDO
	3722	ENDDO
	3723	! ******** MODIF...
	3724
	3725	END SUBROUTINE vertical_diffusion_v_2
	3726
	3727	!=======================================================================
	3728	!=======================================================================
	3729
	3730	SUBROUTINE vertical_diffusion_w_2(tendency, config_flags, mu, &
	3731	defor33, tke, div, xkmv, &
	3732	dn, rdz, &
	3733	ids, ide, jds, jde, kds, kde, &
	3734	ims, ime, jms, jme, kms, kme, &
	3735	its, ite, jts, jte, kts, kte )
	3736
	3737	!-----------------------------------------------------------------------
	3738	! Begin declarations.
	3739
	3740	IMPLICIT NONE
	3741
	3742	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	3743
	3744	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	3745	ims, ime, jms, jme, kms, kme, &
	3746	its, ite, jts, jte, kts, kte
	3747
	3748	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
	3749
	3750	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
	3751
	3752	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
	3753	INTENT(IN ) ::defor33, &
	3754	tke, &
	3755	div, &
	3756	xkmv, &
	3757	rdz
	3758
	3759	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
	3760
	3761	! LOCAL VARS
	3762
	3763	INTEGER :: i, j, k, ktf
	3764
	3765	INTEGER :: i_start, i_end, j_start, j_end
	3766	INTEGER :: is_ext,ie_ext,js_ext,je_ext
	3767
	3768	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
	3769
	3770	! End declarations.
	3771	!-----------------------------------------------------------------------
	3772
	3773	ktf=MIN(kte,kde-1)
	3774
	3775	i_start = its
	3776	i_end = MIN(ite,ide-1)
	3777	j_start = jts
	3778	j_end = MIN(jte,jde-1)
	3779
	3780	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	3781	config_flags%nested) i_start = MAX(ids+1,its)
	3782	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	3783	config_flags%nested) i_end = MIN(ide-2,ite)
	3784	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	3785	config_flags%nested) j_start = MAX(jds+1,jts)
	3786	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	3787	config_flags%nested) j_end = MIN(jde-2,jte)
	3788	IF ( config_flags%periodic_x ) i_start = its
	3789	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
	3790
	3791	! titau3 = titau33
	3792	is_ext=0
	3793	ie_ext=0
	3794	js_ext=0
	3795	je_ext=0
	3796	CALL cal_titau_11_22_33( config_flags, titau3, &
	3797	mu, tke, xkmv, defor33, &
	3798	is_ext, ie_ext, js_ext, je_ext, &
	3799	ids, ide, jds, jde, kds, kde, &
	3800	ims, ime, jms, jme, kms, kme, &
	3801	its, ite, jts, jte, kts, kte )
	3802
	3803	! DO j = j_start, j_end
	3804	! DO k = kts+1, ktf
	3805	! DO i = i_start, i_end
	3806	! titau3(i,k,j)=titau3(i,k,j)*zeta_z(i,j)
	3807	! ENDDO
	3808	! ENDDO
	3809	! ENDDO
	3810
	3811	DO j = j_start, j_end
	3812	DO k = kts+1, ktf
	3813	DO i = i_start, i_end
	3814	tendency(i,k,j)=tendency(i,k,j)-rdz(i,k,j)*(titau3(i,k,j)-titau3(i,k-1,j))
	3815	ENDDO
	3816	ENDDO
	3817	ENDDO
	3818
	3819	END SUBROUTINE vertical_diffusion_w_2
	3820
	3821	!=======================================================================
	3822	!=======================================================================
	3823
	3824	SUBROUTINE vertical_diffusion_s( tendency, config_flags, var, mu, xkhv, &
	3825	dn, dnw, rdz, rdzw, fnm, fnp, &
	3826	doing_tke, &
	3827	ids, ide, jds, jde, kds, kde, &
	3828	ims, ime, jms, jme, kms, kme, &
	3829	its, ite, jts, jte, kts, kte )
	3830
	3831	!-----------------------------------------------------------------------
	3832	! Begin declarations.
	3833
	3834	IMPLICIT NONE
	3835
	3836	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	3837
	3838	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	3839	ims, ime, jms, jme, kms, kme, &
	3840	its, ite, jts, jte, kts, kte
	3841
	3842	LOGICAL, INTENT(IN ) :: doing_tke
	3843
	3844	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
	3845	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
	3846	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
	3847	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
	3848
	3849	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
	3850
	3851	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN) :: xkhv
	3852
	3853	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mu
	3854
	3855	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
	3856	INTENT(IN ) :: var, &
	3857	rdz, &
	3858	rdzw
	3859	! LOCAL VARS
	3860
	3861	INTEGER :: i, j, k, ktf
	3862
	3863	INTEGER :: i_start, i_end, j_start, j_end
	3864
	3865	REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: H3, &
	3866	xkxavg, &
	3867	rravg
	3868
	3869	REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
	3870
	3871	! End declarations.
	3872	!-----------------------------------------------------------------------
	3873
	3874	ktf=MIN(kte,kde-1)
	3875
	3876	i_start = its
	3877	i_end = MIN(ite,ide-1)
	3878	j_start = jts
	3879	j_end = MIN(jte,jde-1)
	3880
	3881	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	3882	config_flags%nested) i_start = MAX(ids+1,its)
	3883	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	3884	config_flags%nested) i_end = MIN(ide-2,ite)
	3885	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	3886	config_flags%nested) j_start = MAX(jds+1,jts)
	3887	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	3888	config_flags%nested) j_end = MIN(jde-2,jte)
	3889	IF ( config_flags%periodic_x ) i_start = its
	3890	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
	3891
	3892	IF (doing_tke) THEN
	3893	DO j = j_start, j_end
	3894	DO k = kts,ktf
	3895	DO i = i_start, i_end
	3896	tmptendf(i,k,j)=tendency(i,k,j)
	3897	ENDDO
	3898	ENDDO
	3899	ENDDO
	3900	ENDIF
	3901
	3902	! H3
	3903
	3904	xkxavg = 0.
	3905
	3906	DO j = j_start, j_end
	3907	DO k = kts+1,ktf
	3908	DO i = i_start, i_end
	3909	xkxavg(i,k,j)=fnm(k)xkhv(i,k,j)+fnp(k)xkhv(i,k-1,j)
	3910	H3(i,k,j)=-xkxavg(i,k,j)(var(i,k,j)-var(i,k-1,j))rdz(i,k,j)
	3911	! H3(i,k,j)=-xkxavg(i,k,j)zeta_z(i,j) &
	3912	! (var(i,k,j)-var(i,k-1,j))/dn(k)
	3913	ENDDO
	3914	ENDDO
	3915	ENDDO
	3916
	3917	DO j = j_start, j_end
	3918	DO i = i_start, i_end
	3919	H3(i,kts,j)=0.
	3920	H3(i,ktf+1,j)=0.
	3921	! H3(i,kts,j)=H3(i,kts+1,j)
	3922	! H3(i,ktf+1,j)=H3(i,ktf,j)
	3923	ENDDO
	3924	ENDDO
	3925
	3926	DO j = j_start, j_end
	3927	DO k = kts,ktf
	3928	DO i = i_start, i_end
	3929	tendency(i,k,j)=tendency(i,k,j) &
	3930	-mu(i,j)(H3(i,k+1,j)-H3(i,k,j))rdzw(i,k,j)
	3931	ENDDO
	3932	ENDDO
	3933	ENDDO
	3934
	3935	IF (doing_tke) THEN
	3936	DO j = j_start, j_end
	3937	DO k = kts,ktf
	3938	DO i = i_start, i_end
	3939	tendency(i,k,j)=tmptendf(i,k,j)+2.* &
	3940	(tendency(i,k,j)-tmptendf(i,k,j))
	3941	ENDDO
	3942	ENDDO
	3943	ENDDO
	3944	ENDIF
	3945
	3946	END SUBROUTINE vertical_diffusion_s
	3947
	3948	!=======================================================================
	3949	!=======================================================================
	3950
	3951	SUBROUTINE cal_titau_11_22_33( config_flags, titau, &
	3952	mu, tke, xkx, defor, &
	3953	is_ext, ie_ext, js_ext, je_ext, &
	3954	ids, ide, jds, jde, kds, kde, &
	3955	ims, ime, jms, jme, kms, kme, &
	3956	its, ite, jts, jte, kts, kte )
	3957
	3958	! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
	3959	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	3960	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
	3961
	3962	! Purpose: This routine calculates stress terms (taus) for use in
	3963	! the calculation of production of TKE by sheared wind
	3964
	3965	! References: Klemp and Wilhelmson (JAS 1978)
	3966	! Deardorff (B-L Meteor 1980)
	3967	! Chen and Dudhia (NCAR WRF physics report 2000)
	3968
	3969	! Key:
	3970
	3971	!-----------------------------------------------------------------------
	3972	! Begin declarations.
	3973
	3974	IMPLICIT NONE
	3975
	3976	TYPE( grid_config_rec_type ), INTENT( IN ) &
	3977	:: config_flags
	3978
	3979	INTEGER, INTENT( IN ) &
	3980	:: ids, ide, jds, jde, kds, kde, &
	3981	ims, ime, jms, jme, kms, kme, &
	3982	its, ite, jts, jte, kts, kte
	3983
	3984	INTEGER, INTENT( IN ) &
	3985	:: is_ext, ie_ext, js_ext, je_ext
	3986
	3987	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
	3988	:: titau
	3989
	3990	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	3991	:: defor, xkx, tke
	3992
	3993	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	3994	:: mu
	3995
	3996	! Local variables.
	3997
	3998	INTEGER &
	3999	:: i, j, k, ktf, i_start, i_end, j_start, j_end
	4000
	4001	! End declarations.
	4002	!-----------------------------------------------------------------------
	4003
	4004	ktf = MIN( kte, kde-1 )
	4005
	4006	i_start = its
	4007	i_end = ite
	4008	j_start = jts
	4009	j_end = jte
	4010
	4011	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	4012	config_flags%nested) i_start = MAX( ids+1, its )
	4013	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	4014	config_flags%nested) i_end = MIN( ide-1, ite )
	4015	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	4016	config_flags%nested) j_start = MAX( jds+1, jts )
	4017	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	4018	config_flags%nested) j_end = MIN( jde-1, jte )
	4019	IF ( config_flags%periodic_x ) i_start = its
	4020	IF ( config_flags%periodic_x ) i_end = ite
	4021
	4022	i_start = i_start - is_ext
	4023	i_end = i_end + ie_ext
	4024	j_start = j_start - js_ext
	4025	j_end = j_end + je_ext
	4026
	4027	IF ( config_flags%km_opt .EQ. 2) THEN
	4028	DO j = j_start,j_end
	4029	DO k = kts,ktf
	4030	DO i = i_start,i_end
	4031	titau(i,k,j) = mu(i,j) * ( - xkx(i,k,j) * ( defor(i,k,j) ) )
	4032	END DO
	4033	END DO
	4034	END DO
	4035	ELSE
	4036	DO j = j_start, j_end
	4037	DO k = kts, ktf
	4038	DO i = i_start, i_end
	4039	titau(i,k,j) = - mu(i,j) * xkx(i,k,j) * defor(i,k,j)
	4040	END DO
	4041	END DO
	4042	END DO
	4043	END IF
	4044
	4045	END SUBROUTINE cal_titau_11_22_33
	4046
	4047	!=======================================================================
	4048	!=======================================================================
	4049
	4050	SUBROUTINE cal_titau_12_21( config_flags, titau, &
	4051	mu, xkx, defor, &
	4052	is_ext, ie_ext, js_ext, je_ext, &
	4053	ids, ide, jds, jde, kds, kde, &
	4054	ims, ime, jms, jme, kms, kme, &
	4055	its, ite, jts, jte, kts, kte )
	4056
	4057	! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR
	4058	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	4059	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
	4060
	4061	! Pusrpose This routine calculates the stress terms (taus) for use in
	4062	! the calculation of production of TKE by sheared wind
	4063
	4064	! References: Klemp and Wilhelmson (JAS 1978)
	4065	! Deardorff (B-L Meteor 1980)
	4066	! Chen and Dudhia (NCAR WRF physics report 2000)
	4067
	4068	! Key:
	4069
	4070	!-----------------------------------------------------------------------
	4071	! Begin declarations.
	4072
	4073	IMPLICIT NONE
	4074
	4075	TYPE( grid_config_rec_type), INTENT( IN ) &
	4076	:: config_flags
	4077
	4078	INTEGER, INTENT( IN ) &
	4079	:: ids, ide, jds, jde, kds, kde, &
	4080	ims, ime, jms, jme, kms, kme, &
	4081	its, ite, jts, jte, kts, kte
	4082
	4083	INTEGER, INTENT( IN ) &
	4084	:: is_ext, ie_ext, js_ext, je_ext
	4085
	4086	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
	4087	:: titau
	4088
	4089	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	4090	:: defor, xkx
	4091
	4092	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	4093	:: mu
	4094
	4095	! Local variables.
	4096
	4097	INTEGER &
	4098	:: i, j, k, ktf, i_start, i_end, j_start, j_end
	4099
	4100	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
	4101	:: xkxavg
	4102
	4103	REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
	4104	:: muavg
	4105
	4106	! End declarations.
	4107	!-----------------------------------------------------------------------
	4108
	4109	ktf = MIN( kte, kde-1 )
	4110
	4111	! Needs one more point in the x and y directions.
	4112
	4113	i_start = its
	4114	i_end = ite
	4115	j_start = jts
	4116	j_end = jte
	4117
	4118	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	4119	config_flags%nested ) i_start = MAX( ids+1, its )
	4120	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	4121	config_flags%nested ) i_end = MIN( ide-1, ite )
	4122	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	4123	config_flags%nested ) j_start = MAX( jds+1, jts )
	4124	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	4125	config_flags%nested ) j_end = MIN( jde-1, jte )
	4126	IF ( config_flags%periodic_x ) i_start = its
	4127	IF ( config_flags%periodic_x ) i_end = ite
	4128
	4129	i_start = i_start - is_ext
	4130	i_end = i_end + ie_ext
	4131	j_start = j_start - js_ext
	4132	j_end = j_end + je_ext
	4133
	4134	DO j = j_start, j_end
	4135	DO k = kts, ktf
	4136	DO i = i_start, i_end
	4137	xkxavg(i,k,j) = 0.25 * ( xkx(i-1,k,j ) + xkx(i,k,j ) + &
	4138	xkx(i-1,k,j-1) + xkx(i,k,j-1) )
	4139	END DO
	4140	END DO
	4141	END DO
	4142
	4143	DO j = j_start, j_end
	4144	DO i = i_start, i_end
	4145	muavg(i,j) = 0.25 * ( mu(i-1,j ) + mu(i,j ) + &
	4146	mu(i-1,j-1) + mu(i,j-1) )
	4147	END DO
	4148	END DO
	4149
	4150	! titau12 or titau21
	4151
	4152	DO j = j_start, j_end
	4153	DO k = kts, ktf
	4154	DO i = i_start, i_end
	4155	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
	4156	END DO
	4157	END DO
	4158	END DO
	4159
	4160	END SUBROUTINE cal_titau_12_21
	4161
	4162	!=======================================================================
	4163
	4164	SUBROUTINE cal_titau_13_31( config_flags, titau, &
	4165	defor, mu, xkx, fnm, fnp, &
	4166	is_ext, ie_ext, js_ext, je_ext, &
	4167	ids, ide, jds, jde, kds, kde, &
	4168	ims, ime, jms, jme, kms, kme, &
	4169	its, ite, jts, jte, kts, kte )
	4170
	4171	! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR
	4172	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	4173	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
	4174
	4175	! Purpose: This routine calculates the stress terms (taus) for use in
	4176	! the calculation of production of TKE by sheared wind
	4177
	4178	! References: Klemp and Wilhelmson (JAS 1978)
	4179	! Deardorff (B-L Meteor 1980)
	4180	! Chen and Dudhia (NCAR WRF physics report 2000)
	4181
	4182	! Key:
	4183
	4184	!-----------------------------------------------------------------------
	4185	! Begin declarations.
	4186
	4187	IMPLICIT NONE
	4188
	4189	TYPE( grid_config_rec_type), INTENT( IN ) &
	4190	:: config_flags
	4191
	4192	INTEGER, INTENT( IN ) &
	4193	:: ids, ide, jds, jde, kds, kde, &
	4194	ims, ime, jms, jme, kms, kme, &
	4195	its, ite, jts, jte, kts, kte
	4196
	4197	INTEGER, INTENT( IN ) &
	4198	:: is_ext, ie_ext, js_ext, je_ext
	4199
	4200	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
	4201	:: fnm, fnp
	4202
	4203	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
	4204	:: titau
	4205
	4206	REAL, DIMENSION( ims:ime, kms:kme, jms:jme), INTENT( IN ) &
	4207	:: defor, xkx
	4208
	4209	REAL, DIMENSION( ims:ime, jms:jme), INTENT( IN ) &
	4210	:: mu
	4211
	4212	! Local variables.
	4213
	4214	INTEGER &
	4215	:: i, j, k, ktf, i_start, i_end, j_start, j_end
	4216
	4217	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
	4218	:: xkxavg
	4219
	4220	REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
	4221	:: muavg
	4222
	4223	! End declarations.
	4224	!-----------------------------------------------------------------------
	4225
	4226	ktf = MIN( kte, kde-1 )
	4227
	4228	! Find ide-1 and jde-1 for averaging to p point.
	4229
	4230	i_start = its
	4231	i_end = ite
	4232	j_start = jts
	4233	j_end = MIN( jte, jde-1 )
	4234
	4235	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	4236	config_flags%nested) i_start = MAX( ids+1, its )
	4237	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	4238	config_flags%nested) i_end = MIN( ide-1, ite )
	4239	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	4240	config_flags%nested) j_start = MAX( jds+1, jts )
	4241	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	4242	config_flags%nested) j_end = MIN( jde-2, jte )
	4243	IF ( config_flags%periodic_x ) i_start = its
	4244	IF ( config_flags%periodic_x ) i_end = ite
	4245
	4246	i_start = i_start - is_ext
	4247	i_end = i_end + ie_ext
	4248	j_start = j_start - js_ext
	4249	j_end = j_end + je_ext
	4250
	4251	DO j = j_start, j_end
	4252	DO k = kts+1, ktf
	4253	DO i = i_start, i_end
	4254	xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i-1,k ,j) ) + &
	4255	fnp(k) * ( xkx(i,k-1,j) + xkx(i-1,k-1,j) ) )
	4256	END DO
	4257	END DO
	4258	END DO
	4259
	4260	DO j = j_start, j_end
	4261	DO i = i_start, i_end
	4262	muavg(i,j) = 0.5 * ( mu(i,j) + mu(i-1,j) )
	4263	END DO
	4264	END DO
	4265
	4266	DO j = j_start, j_end
	4267	DO k = kts+1, ktf
	4268	DO i = i_start, i_end
	4269	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
	4270	ENDDO
	4271	ENDDO
	4272	ENDDO
	4273
	4274	DO j = j_start, j_end
	4275	DO i = i_start, i_end
	4276	titau(i,kts ,j) = 0.0
	4277	titau(i,ktf+1,j) = 0.0
	4278	ENDDO
	4279	ENDDO
	4280
	4281	END SUBROUTINE cal_titau_13_31
	4282
	4283	!=======================================================================
	4284	!=======================================================================
	4285
	4286	SUBROUTINE cal_titau_23_32( config_flags, titau, defor, &
	4287	mu, xkx, fnm, fnp, &
	4288	is_ext, ie_ext, js_ext, je_ext, &
	4289	ids, ide, jds, jde, kds, kde, &
	4290	ims, ime, jms, jme, kms, kme, &
	4291	its, ite, jts, jte, kts, kte )
	4292
	4293	! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
	4294	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	4295	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
	4296
	4297	! Purpose: This routine calculates stress terms (taus) for use in
	4298	! the calculation of production of TKE by sheared wind
	4299
	4300	! References: Klemp and Wilhelmson (JAS 1978)
	4301	! Deardorff (B-L Meteor 1980)
	4302	! Chen and Dudhia (NCAR WRF physics report 2000)
	4303
	4304	! Key:
	4305
	4306	!-----------------------------------------------------------------------
	4307	! Begin declarations.
	4308
	4309	IMPLICIT NONE
	4310
	4311	TYPE( grid_config_rec_type ), INTENT( IN ) &
	4312	:: config_flags
	4313
	4314	INTEGER, INTENT( IN ) &
	4315	:: ids, ide, jds, jde, kds, kde, &
	4316	ims, ime, jms, jme, kms, kme, &
	4317	its, ite, jts, jte, kts, kte
	4318
	4319	INTEGER, INTENT( IN ) &
	4320	:: is_ext,ie_ext,js_ext,je_ext
	4321
	4322	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
	4323	:: fnm, fnp
	4324
	4325	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
	4326	:: titau
	4327
	4328	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	4329	:: defor, xkx
	4330
	4331	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	4332	:: mu
	4333
	4334	! Local variables.
	4335
	4336	INTEGER &
	4337	:: i, j, k, ktf, i_start, i_end, j_start, j_end
	4338
	4339	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
	4340	:: xkxavg
	4341
	4342	REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
	4343	:: muavg
	4344
	4345	! End declarations.
	4346	!-----------------------------------------------------------------------
	4347
	4348	ktf = MIN( kte, kde-1 )
	4349
	4350	! Find ide-1 and jde-1 for averaging to p point.
	4351
	4352	i_start = its
	4353	i_end = MIN( ite, ide-1 )
	4354	j_start = jts
	4355	j_end = jte
	4356
	4357	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	4358	config_flags%nested) i_start = MAX( ids+1, its )
	4359	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	4360	config_flags%nested) i_end = MIN( ide-2, ite )
	4361	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	4362	config_flags%nested) j_start = MAX( jds+1, jts )
	4363	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	4364	config_flags%nested) j_end = MIN( jde-1, jte )
	4365	IF ( config_flags%periodic_x ) i_start = its
	4366	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	4367
	4368	i_start = i_start - is_ext
	4369	i_end = i_end + ie_ext
	4370	j_start = j_start - js_ext
	4371	j_end = j_end + je_ext
	4372
	4373	DO j = j_start, j_end
	4374	DO k = kts+1, ktf
	4375	DO i = i_start, i_end
	4376	xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i,k ,j-1) ) + &
	4377	fnp(k) * ( xkx(i,k-1,j) + xkx(i,k-1,j-1) ) )
	4378	END DO
	4379	END DO
	4380	END DO
	4381
	4382	DO j = j_start, j_end
	4383	DO i = i_start, i_end
	4384	muavg(i,j) = 0.5 * ( mu(i,j) + mu(i,j-1) )
	4385	END DO
	4386	END DO
	4387
	4388	DO j = j_start, j_end
	4389	DO k = kts+1, ktf
	4390	DO i = i_start, i_end
	4391	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
	4392	END DO
	4393	END DO
	4394	END DO
	4395
	4396	DO j = j_start, j_end
	4397	DO i = i_start, i_end
	4398	titau(i,kts ,j) = 0.0
	4399	titau(i,ktf+1,j) = 0.0
	4400	END DO
	4401	END DO
	4402
	4403	END SUBROUTINE cal_titau_23_32
	4404
	4405	!=======================================================================
	4406	!=======================================================================
	4407
	4408	SUBROUTINE phy_bc ( config_flags,div,defor11,defor22,defor33, &
	4409	defor12,defor13,defor23,xkmh,xkmv,xkhh,xkhv,tke, &
	4410	RUBLTEN, RVBLTEN, &
	4411	ids, ide, jds, jde, kds, kde, &
	4412	ims, ime, jms, jme, kms, kme, &
	4413	ips, ipe, jps, jpe, kps, kpe, &
	4414	its, ite, jts, jte, kts, kte )
	4415
	4416	!------------------------------------------------------------------------------
	4417	! Begin declarations.
	4418
	4419	IMPLICIT NONE
	4420
	4421	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
	4422
	4423	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
	4424	ims, ime, jms, jme, kms, kme, &
	4425	ips, ipe, jps, jpe, kps, kpe, &
	4426	its, ite, jts, jte, kts, kte
	4427
	4428	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::RUBLTEN, &
	4429	RVBLTEN, &
	4430	defor11, &
	4431	defor22, &
	4432	defor33, &
	4433	defor12, &
	4434	defor13, &
	4435	defor23, &
	4436	xkmh, &
	4437	xkmv, &
	4438	xkhh, &
	4439	xkhv, &
	4440	tke, &
	4441	div
	4442
	4443	! End declarations.
	4444	!-----------------------------------------------------------------------
	4445
	4446	IF( (config_flags%bl_pbl_physics .GT. 0) &
	4447	.OR. (config_flags%modif_wrf) ) THEN
	4448	!****MARS: Always do that ?? > yes, for periodic simu -- grep in share
	4449
	4450	CALL set_physical_bc3d( RUBLTEN , 't', config_flags, &
	4451	ids, ide, jds, jde, kds, kde, &
	4452	ims, ime, jms, jme, kms, kme, &
	4453	ips, ipe, jps, jpe, kps, kpe, &
	4454	its, ite, jts, jte, kts, kte )
	4455
	4456	CALL set_physical_bc3d( RVBLTEN , 't', config_flags, &
	4457	ids, ide, jds, jde, kds, kde, &
	4458	ims, ime, jms, jme, kms, kme, &
	4459	ips, ipe, jps, jpe, kps, kpe, &
	4460	its, ite, jts, jte, kts, kte )
	4461
	4462	ENDIF
	4463
	4464	! move out of the conditional, below; horiz coeffs needed for
	4465	! all diff_opt cases. JM
	4466
	4467	CALL set_physical_bc3d( xkmh , 't', config_flags, &
	4468	ids, ide, jds, jde, kds, kde, &
	4469	ims, ime, jms, jme, kms, kme, &
	4470	ips, ipe, jps, jpe, kps, kpe, &
	4471	its, ite, jts, jte, kts, kte )
	4472
	4473	CALL set_physical_bc3d( xkhh , 't', config_flags, &
	4474	ids, ide, jds, jde, kds, kde, &
	4475	ims, ime, jms, jme, kms, kme, &
	4476	ips, ipe, jps, jpe, kps, kpe, &
	4477	its, ite, jts, jte, kts, kte )
	4478
	4479	IF(config_flags%diff_opt .eq. 2) THEN
	4480
	4481	CALL set_physical_bc3d( xkmv , 't', config_flags, &
	4482	ids, ide, jds, jde, kds, kde, &
	4483	ims, ime, jms, jme, kms, kme, &
	4484	ips, ipe, jps, jpe, kps, kpe, &
	4485	its, ite, jts, jte, kts, kte )
	4486
	4487	CALL set_physical_bc3d( xkhv , 't', config_flags, &
	4488	ids, ide, jds, jde, kds, kde, &
	4489	ims, ime, jms, jme, kms, kme, &
	4490	ips, ipe, jps, jpe, kps, kpe, &
	4491	its, ite, jts, jte, kts, kte )
	4492
	4493	CALL set_physical_bc3d( tke , 't', config_flags, &
	4494	ids, ide, jds, jde, kds, kde, &
	4495	ims, ime, jms, jme, kms, kme, &
	4496	ips, ipe, jps, jpe, kps, kpe, &
	4497	its, ite, jts, jte, kts, kte )
	4498
	4499	CALL set_physical_bc3d( div , 't', config_flags, &
	4500	ids, ide, jds, jde, kds, kde, &
	4501	ims, ime, jms, jme, kms, kme, &
	4502	ips, ipe, jps, jpe, kps, kpe, &
	4503	its, ite, jts, jte, kts, kte )
	4504
	4505	CALL set_physical_bc3d( defor11 , 't', config_flags, &
	4506	ids, ide, jds, jde, kds, kde, &
	4507	ims, ime, jms, jme, kms, kme, &
	4508	ips, ipe, jps, jpe, kps, kpe, &
	4509	its, ite, jts, jte, kts, kte )
	4510
	4511	CALL set_physical_bc3d( defor22 , 't', config_flags, &
	4512	ids, ide, jds, jde, kds, kde, &
	4513	ims, ime, jms, jme, kms, kme, &
	4514	ips, ipe, jps, jpe, kps, kpe, &
	4515	its, ite, jts, jte, kts, kte )
	4516
	4517	CALL set_physical_bc3d( defor33 , 't', config_flags, &
	4518	ids, ide, jds, jde, kds, kde, &
	4519	ims, ime, jms, jme, kms, kme, &
	4520	ips, ipe, jps, jpe, kps, kpe, &
	4521	its, ite, jts, jte, kts, kte )
	4522
	4523	CALL set_physical_bc3d( defor12 , 'd', config_flags, &
	4524	ids, ide, jds, jde, kds, kde, &
	4525	ims, ime, jms, jme, kms, kme, &
	4526	ips, ipe, jps, jpe, kps, kpe, &
	4527	its, ite, jts, jte, kts, kte )
	4528
	4529	CALL set_physical_bc3d( defor13 , 'e', config_flags, &
	4530	ids, ide, jds, jde, kds, kde, &
	4531	ims, ime, jms, jme, kms, kme, &
	4532	ips, ipe, jps, jpe, kps, kpe, &
	4533	its, ite, jts, jte, kts, kte )
	4534
	4535	CALL set_physical_bc3d( defor23 , 'f', config_flags, &
	4536	ids, ide, jds, jde, kds, kde, &
	4537	ims, ime, jms, jme, kms, kme, &
	4538	ips, ipe, jps, jpe, kps, kpe, &
	4539	its, ite, jts, jte, kts, kte )
	4540
	4541	ENDIF
	4542
	4543	END SUBROUTINE phy_bc
	4544
	4545	!=======================================================================
	4546	!=======================================================================
	4547
	4548	SUBROUTINE tke_rhs( tendency, BN2, config_flags, &
	4549	defor11, defor22, defor33, &
	4550	defor12, defor13, defor23, &
	4551	u, v, w, div, tke, mu, &
	4552	theta, p, p8w, t8w, z, fnm, fnp, &
	4553	cf1, cf2, cf3, msftx, msfty, &
	4554	xkmh, xkmv, xkhv, &
	4555	rdx, rdy, dx, dy, dt, zx, zy, &
	4556	rdz, rdzw, dn, dnw, isotropic, &
	4557	hfx, qfx, qv, ust, rho, &
	4558	ids, ide, jds, jde, kds, kde, &
	4559	ims, ime, jms, jme, kms, kme, &
	4560	its, ite, jts, jte, kts, kte )
	4561
	4562	!-----------------------------------------------------------------------
	4563	! Begin declarations.
	4564
	4565	IMPLICIT NONE
	4566
	4567	TYPE( grid_config_rec_type ), INTENT( IN ) &
	4568	:: config_flags
	4569
	4570	INTEGER, INTENT( IN ) &
	4571	:: ids, ide, jds, jde, kds, kde, &
	4572	ims, ime, jms, jme, kms, kme, &
	4573	its, ite, jts, jte, kts, kte
	4574
	4575	INTEGER, INTENT( IN ) :: isotropic
	4576	REAL, INTENT( IN ) &
	4577	:: cf1, cf2, cf3, dt, rdx, rdy, dx, dy
	4578
	4579	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
	4580	:: fnm, fnp, dnw, dn
	4581
	4582	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	4583	:: msftx, msfty
	4584
	4585	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	4586	:: tendency
	4587
	4588	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	4589	:: defor11, defor22, defor33, defor12, defor13, defor23, &
	4590	div, BN2, tke, xkmh, xkmv, xkhv, zx, zy, u, v, w, theta, &
	4591	p, p8w, t8w, z, rdz, rdzw
	4592
	4593	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	4594	:: mu
	4595
	4596	REAL, DIMENSION ( ims:ime, jms:jme ), INTENT( IN ) &
	4597	:: hfx, ust, qfx
	4598	REAL, DIMENSION ( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) &
	4599	:: qv, rho
	4600
	4601	! Local variables.
	4602
	4603	INTEGER &
	4604	:: i, j, k, ktf, i_start, i_end, j_start, j_end
	4605
	4606	! End declarations.
	4607	!-----------------------------------------------------------------------
	4608
	4609	CALL tke_shear( tendency, config_flags, &
	4610	defor11, defor22, defor33, &
	4611	defor12, defor13, defor23, &
	4612	u, v, w, tke, ust, mu, fnm, fnp, &
	4613	cf1, cf2, cf3, msftx, msfty, &
	4614	xkmh, xkmv, &
	4615	rdx, rdy, zx, zy, rdz, rdzw, dnw, dn, &
	4616	ids, ide, jds, jde, kds, kde, &
	4617	ims, ime, jms, jme, kms, kme, &
	4618	its, ite, jts, jte, kts, kte )
	4619
	4620	CALL tke_buoyancy( tendency, config_flags, mu, &
	4621	tke, xkhv, BN2, theta, dt, &
	4622	hfx, qfx, qv, rho, &
	4623	ids, ide, jds, jde, kds, kde, &
	4624	ims, ime, jms, jme, kms, kme, &
	4625	its, ite, jts, jte, kts, kte )
	4626
	4627	CALL tke_dissip( tendency, config_flags, &
	4628	mu, tke, bn2, theta, p8w, t8w, z, &
	4629	dx, dy,rdz, rdzw, isotropic, &
	4630	msftx, msfty, &
	4631	ids, ide, jds, jde, kds, kde, &
	4632	ims, ime, jms, jme, kms, kme, &
	4633	its, ite, jts, jte, kts, kte )
	4634
	4635	! Set a lower limit on TKE.
	4636
	4637	ktf = MIN( kte, kde-1 )
	4638	i_start = its
	4639	i_end = MIN( ite, ide-1 )
	4640	j_start = jts
	4641	j_end = MIN( jte, jde-1 )
	4642
	4643	IF ( config_flags%open_xs .or. config_flags%specified .or. &
	4644	config_flags%nested) i_start = MAX(ids+1,its)
	4645	IF ( config_flags%open_xe .or. config_flags%specified .or. &
	4646	config_flags%nested) i_end = MIN(ide-2,ite)
	4647	IF ( config_flags%open_ys .or. config_flags%specified .or. &
	4648	config_flags%nested) j_start = MAX(jds+1,jts)
	4649	IF ( config_flags%open_ye .or. config_flags%specified .or. &
	4650	config_flags%nested) j_end = MIN(jde-2,jte)
	4651	IF ( config_flags%periodic_x ) i_start = its
	4652	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	4653
	4654	DO j = j_start, j_end
	4655	DO k = kts, ktf
	4656	DO i = i_start, i_end
	4657	tendency(i,k,j) = max( tendency(i,k,j), -mu(i,j) * max( 0.0 , tke(i,k,j) ) / dt )
	4658	END DO
	4659	END DO
	4660	END DO
	4661
	4662	END SUBROUTINE tke_rhs
	4663
	4664	!=======================================================================
	4665	!=======================================================================
	4666
	4667	SUBROUTINE tke_buoyancy( tendency, config_flags, mu, &
	4668	tke, xkhv, BN2, theta, dt, &
	4669	hfx, qfx, qv, rho, &
	4670	ids, ide, jds, jde, kds, kde, &
	4671	ims, ime, jms, jme, kms, kme, &
	4672	its, ite, jts, jte, kts, kte )
	4673
	4674	!-----------------------------------------------------------------------
	4675	! Begin declarations.
	4676
	4677	IMPLICIT NONE
	4678
	4679	TYPE( grid_config_rec_type ), INTENT( IN ) &
	4680	:: config_flags
	4681
	4682	INTEGER, INTENT( IN ) &
	4683	:: ids, ide, jds, jde, kds, kde, &
	4684	ims, ime, jms, jme, kms, kme, &
	4685	its, ite, jts, jte, kts, kte
	4686
	4687	REAL, INTENT( IN ) &
	4688	:: dt
	4689
	4690	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	4691	:: tendency
	4692
	4693	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	4694	:: xkhv, tke, BN2, theta
	4695
	4696	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	4697	:: mu
	4698
	4699	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) &
	4700	:: qv, rho
	4701
	4702	REAL, DIMENSION(ims:ime, jms:jme ), INTENT ( IN ) :: hfx, qfx
	4703
	4704	! Local variables.
	4705
	4706	INTEGER &
	4707	:: i, j, k, ktf
	4708
	4709	INTEGER &
	4710	:: i_start, i_end, j_start, j_end
	4711
	4712	REAL :: heat_flux, heat_flux0
	4713
	4714	REAL :: cpm
	4715
	4716	! End declarations.
	4717	!-----------------------------------------------------------------------
	4718
	4719	!-----------------------------------------------------------------------
	4720	! Add to the TKE tendency the term that accounts for production of TKE
	4721	! due to buoyant motions.
	4722
	4723	ktf = MIN( kte, kde-1 )
	4724	i_start = its
	4725	i_end = MIN( ite, ide-1 )
	4726	j_start = jts
	4727	j_end = MIN( jte, jde-1 )
	4728
	4729	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	4730	config_flags%nested ) i_start = MAX( ids+1, its )
	4731	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	4732	config_flags%nested ) i_end = MIN( ide-2, ite )
	4733	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	4734	config_flags%nested ) j_start = MAX( jds+1, jts )
	4735	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	4736	config_flags%nested ) j_end = MIN( jde-2, jte )
	4737	IF ( config_flags%periodic_x ) i_start = its
	4738	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	4739
	4740	DO j = j_start, j_end
	4741	DO k = kts+1, ktf
	4742	DO i = i_start, i_end
	4743	tendency(i,k,j) = tendency(i,k,j) - mu(i,j) * xkhv(i,k,j) * BN2(i,k,j)
	4744	END DO
	4745	END DO
	4746	END DO
	4747
	4748	! MARTA: change in the computation of the tke's tendency at the surface.
	4749	! the buoyancy flux is the average of the surface heat flux (0.06) and the
	4750	! flux at the first w level
	4751	!
	4752	! WCS 040331
	4753
	4754	hflux: SELECT CASE( config_flags%isfflx )
	4755	CASE (0,2) ! with fixed surface heat flux given in the namelist
	4756	heat_flux0 = config_flags%tke_heat_flux ! constant heat flux value
	4757	! set in namelist.input
	4758	! LES mods
	4759	K=KTS
	4760	PRINT , 'USE constant H/rho/cp : e.g. value of Hs ', heat_flux0cp*rho(i_start,k,j_start)
	4761	DO j = j_start, j_end
	4762	DO i = i_start, i_end
	4763	heat_flux = heat_flux0
	4764	tendency(i,k,j)= tendency(i,k,j) - &
	4765	mu(i,j)((xkhv(i,k,j)BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.
	4766
	4767	ENDDO
	4768	ENDDO
	4769
	4770	CASE (1) ! use surface heat flux computed from surface routine
	4771	PRINT *, 'USE surface heat flux from PHYSICS', hfx(i_start,j_start)
	4772	K=KTS
	4773	DO j = j_start, j_end
	4774	DO i = i_start, i_end
	4775	cpm = cp * (1. + 0.8*qv(i,k,j))
	4776	heat_flux = (hfx(i,j)/cpm)/rho(i,k,j)
	4777
	4778	tendency(i,k,j)= tendency(i,k,j) - &
	4779	mu(i,j)((xkhv(i,k,j)BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.
	4780
	4781	ENDDO
	4782	ENDDO
	4783
	4784	CASE DEFAULT
	4785	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
	4786	END SELECT hflux
	4787
	4788	! end of MARTA/WCS change
	4789
	4790	! The tendency array now includes production of TKE from buoyant
	4791	! motions.
	4792	!-----------------------------------------------------------------------
	4793
	4794	END SUBROUTINE tke_buoyancy
	4795
	4796	!=======================================================================
	4797	!=======================================================================
	4798
	4799	SUBROUTINE tke_dissip( tendency, config_flags, &
	4800	mu, tke, bn2, theta, p8w, t8w, z, &
	4801	dx, dy, rdz, rdzw, isotropic, &
	4802	msftx, msfty, &
	4803	ids, ide, jds, jde, kds, kde, &
	4804	ims, ime, jms, jme, kms, kme, &
	4805	its, ite, jts, jte, kts, kte )
	4806
	4807	! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
	4808	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	4809	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
	4810
	4811	! Purpose: This routine calculates dissipation of turbulent kinetic
	4812	! energy.
	4813
	4814	! References: Klemp and Wilhelmson (JAS 1978)
	4815	! Deardorff (B-L Meteor 1980)
	4816	! Chen and Dudhia (NCAR WRF physics report 2000)
	4817
	4818	!-----------------------------------------------------------------------
	4819	! Begin declarations.
	4820
	4821	IMPLICIT NONE
	4822
	4823	TYPE( grid_config_rec_type ), INTENT( IN ) &
	4824	:: config_flags
	4825
	4826	INTEGER, INTENT( IN ) &
	4827	:: ids, ide, jds, jde, kds, kde, &
	4828	ims, ime, jms, jme, kms, kme, &
	4829	its, ite, jts, jte, kts, kte
	4830
	4831	INTEGER, INTENT( IN ) :: isotropic
	4832	REAL, INTENT( IN ) &
	4833	:: dx, dy
	4834
	4835	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	4836	:: tendency
	4837
	4838	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	4839	:: tke, bn2, theta, p8w, t8w, z, rdz, rdzw
	4840
	4841	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	4842	:: mu
	4843
	4844	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	4845	:: msftx, msfty
	4846	! Local variables.
	4847
	4848	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
	4849	:: dthrdn
	4850
	4851	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
	4852	:: l_scale
	4853
	4854	REAL, DIMENSION( its:ite ) &
	4855	:: sumtke, sumtkez
	4856
	4857	INTEGER &
	4858	:: i, j, k, ktf, i_start, i_end, j_start, j_end
	4859
	4860	REAL &
	4861	:: disp_len, deltas, coefc, tmpdz, len_s, thetasfc, &
	4862	thetatop, len_0, tketmp, tmp, ce1, ce2, c_k
	4863
	4864	! End declarations.
	4865	!-----------------------------------------------------------------------
	4866	c_k = config_flags%c_k
	4867
	4868	ce1 = ( c_k / 0.10 ) * 0.19
	4869	ce2 = max( 0.0 , 0.93 - ce1 )
	4870
	4871	ktf = MIN( kte, kde-1 )
	4872	i_start = its
	4873	i_end = MIN(ite,ide-1)
	4874	j_start = jts
	4875	j_end = MIN(jte,jde-1)
	4876
	4877	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	4878	config_flags%nested) i_start = MAX( ids+1, its )
	4879	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	4880	config_flags%nested) i_end = MIN( ide-2, ite )
	4881	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	4882	config_flags%nested) j_start = MAX( jds+1, jts )
	4883	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	4884	config_flags%nested) j_end = MIN( jde-2, jte )
	4885	IF ( config_flags%periodic_x ) i_start = its
	4886	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	4887
	4888	CALL calc_l_scale( config_flags, tke, BN2, l_scale, &
	4889	i_start, i_end, ktf, j_start, j_end, &
	4890	dx, dy, rdzw, msftx, msfty, &
	4891	ids, ide, jds, jde, kds, kde, &
	4892	ims, ime, jms, jme, kms, kme, &
	4893	its, ite, jts, jte, kts, kte )
	4894	DO j = j_start, j_end
	4895	DO k = kts, ktf
	4896	DO i = i_start, i_end
	4897	deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
	4898	tketmp = MAX( tke(i,k,j), 1.0e-6 )
	4899
	4900	! Apply Deardorff's (1980) "wall effect" at the bottom of the domain.
	4901	! For LES with fine grid, no need for this wall effect!
	4902
	4903	IF ( k .eq. kts .or. k .eq. ktf ) then
	4904	coefc = 3.9
	4905	ELSE
	4906	coefc = ce1 + ce2 * l_scale(i,k,j) / deltas
	4907	END IF
	4908
	4909	tendency(i,k,j) = tendency(i,k,j) - &
	4910	mu(i,j) * coefc * tketmp**1.5 / l_scale(i,k,j)
	4911	END DO
	4912	END DO
	4913	END DO
	4914
	4915	END SUBROUTINE tke_dissip
	4916
	4917	!=======================================================================
	4918	!=======================================================================
	4919
	4920	SUBROUTINE tke_shear( tendency, config_flags, &
	4921	defor11, defor22, defor33, &
	4922	defor12, defor13, defor23, &
	4923	u, v, w, tke, ust, mu, fnm, fnp, &
	4924	cf1, cf2, cf3, msftx, msfty, &
	4925	xkmh, xkmv, &
	4926	rdx, rdy, zx, zy, rdz, rdzw, dn, dnw, &
	4927	ids, ide, jds, jde, kds, kde, &
	4928	ims, ime, jms, jme, kms, kme, &
	4929	its, ite, jts, jte, kts, kte )
	4930
	4931	! History: Sep 2003 Rewritten by George Bryan and Jason Knievel,
	4932	! NCAR
	4933	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
	4934	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
	4935
	4936	! Purpose: This routine calculates the production of turbulent
	4937	! kinetic energy by stresses due to sheared wind.
	4938
	4939	! References: Klemp and Wilhelmson (JAS 1978)
	4940	! Deardorff (B-L Meteor 1980)
	4941	! Chen and Dudhia (NCAR WRF physics report 2000)
	4942
	4943	! Key:
	4944
	4945	! avg temporary working array
	4946	! cf1
	4947	! cf2
	4948	! cf3
	4949	! defor11 deformation term ( du/dx + du/dx )
	4950	! defor12 deformation term ( dv/dx + du/dy ); same as defor21
	4951	! defor13 deformation term ( dw/dx + du/dz ); same as defor31
	4952	! defor22 deformation term ( dv/dy + dv/dy )
	4953	! defor23 deformation term ( dw/dy + dv/dz ); same as defor32
	4954	! defor33 deformation term ( dw/dz + dw/dz )
	4955	! div 3-d divergence
	4956	! dn
	4957	! dnw
	4958	! fnm
	4959	! fnp
	4960	! msftx
	4961	! msfty
	4962	! rdx
	4963	! rdy
	4964	! tendency
	4965	! titau tau (stress tensor) with a tilde, indicating division by
	4966	! a map-scale factor and the fraction of the total modeled
	4967	! atmosphere beneath a given altitude (titau = tau/m/zeta)
	4968	! tke turbulent kinetic energy
	4969
	4970	!-----------------------------------------------------------------------
	4971	! Begin declarations.
	4972
	4973	IMPLICIT NONE
	4974
	4975	TYPE( grid_config_rec_type ), INTENT( IN ) &
	4976	:: config_flags
	4977
	4978	INTEGER, INTENT( IN ) &
	4979	:: ids, ide, jds, jde, kds, kde, &
	4980	ims, ime, jms, jme, kms, kme, &
	4981	its, ite, jts, jte, kts, kte
	4982
	4983	REAL, INTENT( IN ) &
	4984	:: cf1, cf2, cf3, rdx, rdy
	4985
	4986	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
	4987	:: fnm, fnp, dn, dnw
	4988
	4989	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	4990	:: msftx, msfty
	4991
	4992	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
	4993	:: tendency
	4994
	4995	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	4996	:: defor11, defor22, defor33, defor12, defor13, defor23, &
	4997	tke, xkmh, xkmv, zx, zy, u, v, w, rdz, rdzw
	4998
	4999	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	5000	:: mu
	5001
	5002	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
	5003	:: ust
	5004
	5005	! Local variables.
	5006
	5007	INTEGER &
	5008	:: i, j, k, ktf, ktes1, ktes2, &
	5009	i_start, i_end, j_start, j_end, &
	5010	is_ext, ie_ext, js_ext, je_ext
	5011
	5012	REAL &
	5013	:: mtau
	5014
	5015	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
	5016	:: avg, titau, tmp2
	5017
	5018	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
	5019	:: titau12, tmp1, zxavg, zyavg
	5020
	5021	REAL :: absU, cd0, Cd
	5022
	5023	! End declarations.
	5024	!-----------------------------------------------------------------------
	5025
	5026	ktf = MIN( kte, kde-1 )
	5027	ktes1 = kte-1
	5028	ktes2 = kte-2
	5029
	5030	i_start = its
	5031	i_end = MIN( ite, ide-1 )
	5032	j_start = jts
	5033	j_end = MIN( jte, jde-1 )
	5034
	5035	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
	5036	config_flags%nested ) i_start = MAX( ids+1, its )
	5037	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
	5038	config_flags%nested ) i_end = MIN( ide-2, ite )
	5039	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
	5040	config_flags%nested ) j_start = MAX( jds+1, jts )
	5041	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
	5042	config_flags%nested ) j_end = MIN( jde-2, jte )
	5043	IF ( config_flags%periodic_x ) i_start = its
	5044	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
	5045
	5046	DO j = j_start, j_end
	5047	DO k = kts, ktf
	5048	DO i = i_start, i_end
	5049	zxavg(i,k,j) = 0.25 * ( zx(i,k ,j) + zx(i+1,k ,j) + &
	5050	zx(i,k+1,j) + zx(i+1,k+1,j) )
	5051	zyavg(i,k,j) = 0.25 * ( zy(i,k ,j) + zy(i,k ,j+1) + &
	5052	zy(i,k+1,j) + zy(i,k+1,j+1) )
	5053	END DO
	5054	END DO
	5055	END DO
	5056
	5057	! Begin calculating production of turbulence due to shear. The approach
	5058	! is to add together contributions from six terms, each of which is the
	5059	! square of a deformation that is then multiplied by an exchange
	5060	! coefficiant. The same exchange coefficient is assumed for horizontal
	5061	! and vertical coefficients for some of the terms (the vertical value is
	5062	! the one used).
	5063
	5064	! For defor11.
	5065
	5066	DO j = j_start, j_end
	5067	DO k = kts, ktf
	5068	DO i = i_start, i_end
	5069	tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
	5070	mu(i,j) * xkmh(i,k,j) * ( ( defor11(i,k,j) )**2 )
	5071	END DO
	5072	END DO
	5073	END DO
	5074
	5075	! For defor22.
	5076
	5077	DO j = j_start, j_end
	5078	DO k = kts, ktf
	5079	DO i = i_start, i_end
	5080	tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
	5081	mu(i,j) * xkmh(i,k,j) * ( ( defor22(i,k,j) )**2 )
	5082	END DO
	5083	END DO
	5084	END DO
	5085
	5086	! For defor33.
	5087
	5088	DO j = j_start, j_end
	5089	DO k = kts, ktf
	5090	DO i = i_start, i_end
	5091	tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
	5092	mu(i,j) * xkmv(i,k,j) * ( ( defor33(i,k,j) )**2 )
	5093	END DO
	5094	END DO
	5095	END DO
	5096
	5097	! For defor12.
	5098
	5099	DO j = j_start, j_end
	5100	DO k = kts, ktf
	5101	DO i = i_start, i_end
	5102	avg(i,k,j) = 0.25 * &
	5103	( ( defor12(i ,k,j)2 ) + ( defor12(i ,k,j+1)2 ) + &
	5104	( defor12(i+1,k,j)2 ) + ( defor12(i+1,k,j+1)2 ) )
	5105	END DO
	5106	END DO
	5107	END DO
	5108
	5109	DO j = j_start, j_end
	5110	DO k = kts, ktf
	5111	DO i = i_start, i_end
	5112	tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmh(i,k,j) * avg(i,k,j)
	5113	END DO
	5114	END DO
	5115	END DO
	5116
	5117	! For defor13.
	5118
	5119	DO j = j_start, j_end
	5120	DO k = kts+1, ktf
	5121	DO i = i_start, i_end+1
	5122	tmp2(i,k,j) = defor13(i,k,j)
	5123	END DO
	5124	END DO
	5125	END DO
	5126
	5127	DO j = j_start, j_end
	5128	DO i = i_start, i_end+1
	5129	tmp2(i,kts ,j) = 0.0
	5130	tmp2(i,ktf+1,j) = 0.0
	5131	END DO
	5132	END DO
	5133
	5134	DO j = j_start, j_end
	5135	DO k = kts, ktf
	5136	DO i = i_start, i_end
	5137	avg(i,k,j) = 0.25 * &
	5138	( ( tmp2(i ,k+1,j)2 ) + ( tmp2(i ,k,j)2 ) + &
	5139	( tmp2(i+1,k+1,j)2 ) + ( tmp2(i+1,k,j)2 ) )
	5140	END DO
	5141	END DO
	5142	END DO
	5143
	5144	DO j = j_start, j_end
	5145	DO k = kts, ktf
	5146	DO i = i_start, i_end
	5147	tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
	5148	END DO
	5149	END DO
	5150	END DO
	5151
	5152	!MARTA: add the drag at the surface; WCS 040331
	5153	K=KTS
	5154
	5155	uflux: SELECT CASE( config_flags%isfflx )
	5156	CASE (0) ! Assume cd a constant, specified in namelist
	5157
	5158	cd0 = config_flags%tke_drag_coefficient ! drag coefficient set
	5159	! in namelist.input
	5160	DO j = j_start, j_end
	5161	DO i = i_start, i_end
	5162
	5163	absU=0.5sqrt((u(i,k,j)+u(i+1,k,j))2+(v(i,k,j)+v(i,k,j+1))*2)
	5164	Cd = cd0
	5165	tendency(i,k,j) = tendency(i,k,j) + &
	5166	mu(i,j)( (u(i,k,j)+u(i+1,k,j))0.5* &
	5167	CdabsU(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 )
	5168
	5169	END DO
	5170	END DO
	5171
	5172	CASE (1,2) ! ustar computed from surface routine
	5173
	5174	PRINT *, 'USE u star (surf drag) from PHYSICS', ust(i_start,j_start)
	5175	DO j = j_start, j_end
	5176	DO i = i_start, i_end
	5177
	5178	absU=0.5sqrt((u(i,k,j)+u(i+1,k,j))2+(v(i,k,j)+v(i,k,j+1))*2)+epsilon
	5179	Cd = (ust(i,j)2)/(absU2)
	5180	tendency(i,k,j) = tendency(i,k,j) + &
	5181	mu(i,j)( (u(i,k,j)+u(i+1,k,j))0.5* &
	5182	CdabsU(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 )
	5183
	5184	END DO
	5185	END DO
	5186
	5187	CASE DEFAULT
	5188	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
	5189	END SELECT uflux
	5190	! end of MARTA/WCS change
	5191
	5192	! For defor23.
	5193
	5194	DO j = j_start, j_end+1
	5195	DO k = kts+1, ktf
	5196	DO i = i_start, i_end
	5197	tmp2(i,k,j) = defor23(i,k,j)
	5198	END DO
	5199	END DO
	5200	END DO
	5201
	5202	DO j = j_start, j_end+1
	5203	DO i = i_start, i_end
	5204	tmp2(i,kts, j) = 0.0
	5205	tmp2(i,ktf+1,j) = 0.0
	5206	END DO
	5207	END DO
	5208
	5209	DO j = j_start, j_end
	5210	DO k = kts, ktf
	5211	DO i = i_start, i_end
	5212	avg(i,k,j) = 0.25 * &
	5213	( ( tmp2(i,k+1,j )2 ) + ( tmp2(i,k,j )2) + &
	5214	( tmp2(i,k+1,j+1)2 ) + ( tmp2(i,k,j+1)2) )
	5215	END DO
	5216	END DO
	5217	END DO
	5218
	5219	DO j = j_start, j_end
	5220	DO k = kts, ktf
	5221	DO i = i_start, i_end
	5222	tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
	5223	END DO
	5224	END DO
	5225	END DO
	5226
	5227	!MARTA: add the drag at the surface; WCS 040331
	5228	K=KTS
	5229
	5230	vflux: SELECT CASE( config_flags%isfflx )
	5231	CASE (0) ! Assume cd a constant, specified in namelist
	5232
	5233	cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient
	5234	! set in namelist.input
	5235	DO j = j_start, j_end
	5236	DO i = i_start, i_end
	5237
	5238	absU=0.5sqrt((u(i,k,j)+u(i+1,k,j))2+(v(i,k,j)+v(i,k,j+1))*2)
	5239	Cd = cd0
	5240	tendency(i,k,j) = tendency(i,k,j) + &
	5241	mu(i,j)( (v(i,k,j)+v(i,k,j+1))0.5* &
	5242	CdabsU(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 )
	5243
	5244	END DO
	5245	END DO
	5246
	5247	CASE (1,2) ! ustar computed from surface routine
	5248
	5249	DO j = j_start, j_end
	5250	DO i = i_start, i_end
	5251
	5252	absU=0.5sqrt((u(i,k,j)+u(i+1,k,j))2+(v(i,k,j)+v(i,k,j+1))*2)+epsilon
	5253	Cd = (ust(i,j)2)/(absU2)
	5254	tendency(i,k,j) = tendency(i,k,j) + &
	5255	mu(i,j)( (v(i,k,j)+v(i,k,j+1))0.5* &
	5256	CdabsU(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 )
	5257
	5258	END DO
	5259	END DO
	5260
	5261	CASE DEFAULT
	5262	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
	5263	END SELECT vflux
	5264	! end of MARTA/WCS change
	5265
	5266	END SUBROUTINE tke_shear
	5267
	5268	!=======================================================================
	5269	!=======================================================================
	5270
	5271	SUBROUTINE compute_diff_metrics( config_flags, ph, phb, z, rdz, rdzw, &
	5272	zx, zy, rdx, rdy, &
	5273	ids, ide, jds, jde, kds, kde, &
	5274	ims, ime, jms, jme, kms, kme, &
	5275	its, ite, jts, jte, kts, kte )
	5276
	5277	!-----------------------------------------------------------------------
	5278	! Begin declarations.
	5279
	5280	IMPLICIT NONE
	5281
	5282	TYPE( grid_config_rec_type ), INTENT( IN ) &
	5283	:: config_flags
	5284
	5285	INTEGER, INTENT( IN ) &
	5286	:: ids, ide, jds, jde, kds, kde, &
	5287	ims, ime, jms, jme, kms, kme, &
	5288	its, ite, jts, jte, kts, kte
	5289
	5290	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
	5291	:: ph, phb
	5292
	5293	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
	5294	:: rdz, rdzw, zx, zy, z
	5295
	5296	REAL, INTENT( IN ) &
	5297	:: rdx, rdy
	5298
	5299	! Local variables.
	5300
	5301	INTEGER &
	5302	:: i, j, k, i_start, i_end, j_start, j_end, ktf
	5303
	5304	! End declarations.
	5305	!-----------------------------------------------------------------------
	5306
	5307	ktf = MIN( kte, kde-1 )
	5308
	5309	! Bug fix, WCS, 22 april 2002.
	5310
	5311	! We need rdzw in halo for average to u and v points.
	5312
	5313	j_start = jts-1
	5314	j_end = jte
	5315
	5316	! Begin with dz computations.
	5317
	5318	DO j = j_start, j_end
	5319
	5320	IF ( ( j_start >= jts ) .AND. ( j_end <= MIN( jte, jde-1 ) ) ) THEN
	5321	i_start = its-1
	5322	i_end = ite
	5323	ELSE
	5324	i_start = its
	5325	i_end = MIN( ite, ide-1 )
	5326	END IF
	5327
	5328	! Compute z at w points for rdz and rdzw computations. We'll switch z
	5329	! to z at p points before returning
	5330
	5331	DO k = 1, kte
	5332
	5333	! Bug fix, WCS, 22 april 2002
	5334
	5335	DO i = i_start, i_end
	5336	z(i,k,j) = ( ph(i,k,j) + phb(i,k,j) ) / g
	5337	END DO
	5338	END DO
	5339
	5340	DO k = 1, ktf
	5341	DO i = i_start, i_end
	5342	rdzw(i,k,j) = 1.0 / ( z(i,k+1,j) - z(i,k,j) )
	5343	END DO
	5344	END DO
	5345
	5346	DO k = 2, ktf
	5347	DO i = i_start, i_end
	5348	rdz(i,k,j) = 2.0 / ( z(i,k+1,j) - z(i,k-1,j) )
	5349	END DO
	5350	END DO
	5351
	5352	! Bug fix, WCS, 22 april 2002; added the following code
	5353
	5354	DO i = i_start, i_end
	5355	rdz(i,1,j) = 2./(z(i,2,j)-z(i,1,j))
	5356	END DO
	5357
	5358	END DO
	5359
	5360	! End bug fix.
	5361
	5362	! Now compute zx and zy; we'll assume that the halo for ph and phb is
	5363	! properly filled.
	5364
	5365	i_start = its
	5366	i_end = MIN( ite, ide-1 )
	5367	j_start = jts
	5368	j_end = MIN( jte, jde-1 )
	5369
	5370	DO j = j_start, j_end
	5371	DO k = 1, kte
	5372	DO i = MAX( ids+1, its ), i_end
	5373	zx(i,k,j) = rdx * ( phb(i,k,j) - phb(i-1,k,j) ) / g
	5374	END DO
	5375	END DO
	5376	END DO
	5377
	5378	DO j = j_start, j_end
	5379	DO k = 1, kte
	5380	DO i = MAX( ids+1, its ), i_end
	5381	zx(i,k,j) = zx(i,k,j) + rdx * ( ph(i,k,j) - ph(i-1,k,j) ) / g
	5382	END DO
	5383	END DO
	5384	END DO
	5385
	5386	DO j = MAX( jds+1, jts ), j_end
	5387	DO k = 1, kte
	5388	DO i = i_start, i_end
	5389	zy(i,k,j) = rdy * ( phb(i,k,j) - phb(i,k,j-1) ) / g
	5390	END DO
	5391	END DO
	5392	END DO
	5393
	5394	DO j = MAX( jds+1, jts ), j_end
	5395	DO k = 1, kte
	5396	DO i = i_start, i_end
	5397	zy(i,k,j) = zy(i,k,j) + rdy * ( ph(i,k,j) - ph(i,k,j-1) ) / g
	5398	END DO
	5399	END DO
	5400	END DO
	5401
	5402	! Some b.c. on zx and zy.
	5403
	5404	IF ( .NOT. config_flags%periodic_x ) THEN
	5405
	5406	IF ( ite == ide ) THEN
	5407	DO j = j_start, j_end
	5408	DO k = 1, ktf
	5409	zx(ide,k,j) = 0.0
	5410	END DO
	5411	END DO
	5412	END IF
	5413
	5414	IF ( its == ids ) THEN
	5415	DO j = j_start, j_end
	5416	DO k = 1, ktf
	5417	zx(ids,k,j) = 0.0
	5418	END DO
	5419	END DO
	5420	END IF
	5421
	5422	ELSE
	5423
	5424	IF ( ite == ide ) THEN
	5425	DO j=j_start,j_end
	5426	DO k=1,ktf
	5427	zx(ide,k,j) = rdx * ( phb(ide,k,j) - phb(ide-1,k,j) ) / g
	5428	END DO
	5429	END DO
	5430
	5431	DO j = j_start, j_end
	5432	DO k = 1, ktf
	5433	zx(ide,k,j) = zx(ide,k,j) + rdx * ( ph(ide,k,j) - ph(ide-1,k,j) ) / g
	5434	END DO
	5435	END DO
	5436	END IF
	5437
	5438	IF ( its == ids ) THEN
	5439	DO j = j_start, j_end
	5440	DO k = 1, ktf
	5441	zx(ids,k,j) = rdx * ( phb(ids,k,j) - phb(ids-1,k,j) ) / g
	5442	END DO
	5443	END DO
	5444
	5445	DO j =j_start,j_end
	5446	DO k =1,ktf
	5447	zx(ids,k,j) = zx(ids,k,j) + rdx * ( ph(ids,k,j) - ph(ids-1,k,j) ) / g
	5448	END DO
	5449	END DO
	5450	END IF
	5451
	5452	END IF
	5453
	5454	IF ( .NOT. config_flags%periodic_y ) THEN
	5455
	5456	IF ( jte == jde ) THEN
	5457	DO k =1, ktf
	5458	DO i =i_start, i_end
	5459	zy(i,k,jde) = 0.0
	5460	END DO
	5461	END DO
	5462	END IF
	5463
	5464	IF ( jts == jds ) THEN
	5465	DO k =1, ktf
	5466	DO i =i_start, i_end
	5467	zy(i,k,jds) = 0.0
	5468	END DO
	5469	END DO
	5470	END IF
	5471
	5472	ELSE
	5473
	5474	IF ( jte == jde ) THEN
	5475	DO j=j_start, j_end
	5476	DO k=1, ktf
	5477	zy(i,k,jde) = rdy * ( phb(i,k,jde) - phb(i,k,jde-1) ) / g
	5478	END DO
	5479	END DO
	5480
	5481	DO j = j_start, j_end
	5482	DO k = 1, ktf
	5483	zy(i,k,jde) = zy(i,k,jde) + rdy * ( ph(i,k,jde) - ph(i,k,jde-1) ) / g
	5484	END DO
	5485	END DO
	5486	END IF
	5487
	5488	IF ( jts == jds ) THEN
	5489	DO j = j_start, j_end
	5490	DO k = 1, ktf
	5491	zy(i,k,jds) = rdy * ( phb(i,k,jds) - phb(i,k,jds-1) ) / g
	5492	END DO
	5493	END DO
	5494
	5495	DO j = j_start, j_end
	5496	DO k = 1, ktf
	5497	zy(i,k,jds) = zy(i,k,jds) + rdy * ( ph(i,k,jds) - ph(i,k,jds-1) ) / g
	5498	END DO
	5499	END DO
	5500	END IF
	5501
	5502	END IF
	5503
	5504	! Calculate z at p points.
	5505
	5506	DO j = j_start, j_end
	5507	DO k = 1, ktf
	5508	DO i = i_start, i_end
	5509	z(i,k,j) = 0.5 * &
	5510	( ph(i,k,j) + phb(i,k,j) + ph(i,k+1,j) + phb(i,k+1,j) ) / g
	5511	END DO
	5512	END DO
	5513	END DO
	5514
	5515	END SUBROUTINE compute_diff_metrics
	5516
	5517	!=======================================================================
	5518	!=======================================================================
	5519
	5520	END MODULE module_diffusion_em
	5521
	5522	!=======================================================================
	5523	!=======================================================================
	5524
	5525

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

Download in other formats:

Original Format