Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_bl_acm.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
Property svn:executable set to ``*
File size: 43.4 KB

Rev	Line
[2759]	1	!WRF:MODEL_LAYER:PHYSICS
	2	!
	3	MODULE module_bl_acm
	4
	5	! USE module_model_constants
	6
	7	REAL, PARAMETER :: RIC = 0.25 ! critical Richardson number
	8	REAL, PARAMETER :: CRANKP = 0.5 ! CRANK-NIC PARAMETER
	9
	10	CONTAINS
	11
	12	!-------------------------------------------------------------------
	13	SUBROUTINE ACMPBL(XTIME, DTPBL, ZNW, SIGMAH, &
	14	U3D, V3D, PP3D, DZ8W, TH3D, T3D, &
	15	QV3D, QC3D, QI3D, RR3D, & !! Moisture
	16	UST, HFX, QFX, TSK, &
	17	PSFC, EP1, G, &
	18	ROVCP, RD, CPD, &
	19	PBLH, KPBL2D, REGIME, &
	20	GZ1OZ0, WSPD, PSIM, MUT, &
	21	RUBLTEN, RVBLTEN, RTHBLTEN, &
	22	RQVBLTEN, RQCBLTEN, RQIBLTEN, &
	23	ids,ide, jds,jde, kds,kde, &
	24	ims,ime, jms,jme, kms,kme, &
	25	its,ite, jts,jte, kts,kte)
	26	!**********************************************************************
	27
	28	! THIS MODULE COMPUTES VERTICAL MIXING IN AND ABOVE THE PBL ACCORDING TO
	29	! THE ASYMMETRICAL CONVECTIVE MODEL, VERSION 2 (ACM2), WHICH IS A COMBINED
	30	! LOCAL NON-LOCAL CLOSURE SCHEME BASED ON THE ORIGINAL ACM (PLEIM AND CHANG 1992)
	31	!
	32	! REFERENCES:
	33	! Pleim (2007) A combined local and non-local closure model for the atmospheric
	34	! boundary layer. Part1: Model description and testing.
	35	! JAMC, 46, 1383-1395
	36	! Pleim (2007) A combined local and non-local closure model for the atmospheric
	37	! boundary layer. Part2: Application and evaluation in a mesoscale
	38	! meteorology model. JAMC, 46, 1396-1409
	39	!
	40	! REVISION HISTORY:
	41	! AX 3/2005 - developed WRF version based on the MM5 PX LSM
	42	! RG and JP 7/2006 - Finished WRF adaptation
	43	!
	44	!**********************************************************************
	45	! ARGUMENT LIST:
	46	!
	47	!... Inputs:
	48	!-- XTIME Time since simulation start (min)
	49	!-- DTPBL PBL time step
	50	!-- ZNW Sigma at full layer
	51	!-- SIGMAH Sigma at half layer
	52	!-- U3D 3D u-velocity interpolated to theta points (m/s)
	53	!-- V3D 3D v-velocity interpolated to theta points (m/s)
	54	!-- PP3D Pressure at half levels (Pa)
	55	!-- DZ8W dz between full levels (m)
	56	!-- TH3D Potential Temperature (K)
	57	!-- T3D Temperature (K)
	58	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
	59	!-- QC3D 3D cloud mixing ratio (Kg/Kg)
	60	!-- QI3D 3D ice mixing ratio (Kg/Kg)
	61	!-- RR3D 3D dry air density (kg/m^3)
	62	!-- UST Friction Velocity (m/s)
	63	!-- HFX Upward heat flux at the surface (w/m^2)
	64	!-- QFX Upward moisture flux at the surface (Kg/m^2/s)
	65	!-- TSK Surface temperature (K)
	66	!-- PSFC Pressure at the surface (Pa)
	67	!-- EP1 Constant for virtual temperature (r_v/r_d-1) (dimensionless)
	68	!-- G Gravity (m/s^2)
	69	!-- ROVCP r/cp
	70	!-- RD gas constant for dry air (j/kg/k)
	71	!-- CPD heat capacity at constant pressure for dry air (j/kg/k)
	72	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
	73	!-- WSPD wind speed at lowest model level (m/s)
	74	!-- PSIM similarity stability function for momentum
	75	!-- MUT Total Mu : Psfc - Ptop
	76
	77	!-- ids start index for i in domain
	78	!-- ide end index for i in domain
	79	!-- jds start index for j in domain
	80	!-- jde end index for j in domain
	81	!-- kds start index for k in domain
	82	!-- kde end index for k in domain
	83	!-- ims start index for i in memory
	84	!-- ime end index for i in memory
	85	!-- jms start index for j in memory
	86	!-- jme end index for j in memory
	87	!-- kms start index for k in memory
	88	!-- kme end index for k in memory
	89	!-- jts start index for j in tile
	90	!-- jte end index for j in tile
	91	!-- kts start index for k in tile
	92	!-- kte end index for k in tile
	93	!
	94	!... Outputs:
	95	!-- PBLH PBL height (m)
	96	!-- KPBL2D K index for PBL layer
	97	!-- REGIME Flag indicating PBL regime (stable, unstable, etc.)
	98	!-- RUBLTEN U tendency due to PBL parameterization (m/s^2)
	99	!-- RVBLTEN V tendency due to PBL parameterization (m/s^2)
	100	!-- RTHBLTEN Theta tendency due to PBL parameterization (K/s)
	101	!-- RQVBLTEN Qv tendency due to PBL parameterization (kg/kg/s)
	102	!-- RQCBLTEN Qc tendency due to PBL parameterization (kg/kg/s)
	103	!-- RQIBLTEN Qi tendency due to PBL parameterization (kg/kg/s)
	104	!
	105	!-------------------------------------------------------------------
	106	IMPLICIT NONE
	107
	108	!.......Arguments
	109	! DECLARATIONS - INTEGER
	110	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	111	ims,ime, jms,jme, kms,kme, &
	112	its,ite, jts,jte, kts,kte, XTIME
	113
	114	! DECLARATIONS - REAL
	115	REAL, INTENT(IN) :: DTPBL, EP1, &
	116	G, ROVCP, RD, CPD
	117	REAL, DIMENSION( kms:kme ), INTENT(IN) :: ZNW, SIGMAH
	118	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
	119	INTENT(IN) :: U3D, V3D, &
	120	PP3D, DZ8W, T3D, &
	121	QV3D, QC3D, QI3D, &
	122	RR3D, TH3D
	123	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST
	124	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
	125	HFX, QFX, TSK, &
	126	PSFC
	127	INTEGER, DIMENSION( ims:ime, jms:jme ), &
	128	INTENT(OUT ) :: KPBL2D
	129	REAL, DIMENSION( ims:ime, jms:jme ), &
	130	INTENT(INOUT) :: PBLH
	131	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: REGIME
	132	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
	133	psim, &
	134	gz1oz0, &
	135	wspd , mut
	136	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
	137	INTENT(INOUT) :: RUBLTEN, RVBLTEN, &
	138	RTHBLTEN, RQVBLTEN, &
	139	RQCBLTEN, RQIBLTEN
	140
	141	!... Local Variables
	142
	143	!... Integer
	144	INTEGER :: J, K
	145	!... Real
	146	REAL, DIMENSION( kts:kte ) :: DSIGH, DSIGHI, DSIGFI
	147	REAL, DIMENSION( 0:kte ) :: SIGMAF
	148	REAL RDT
	149	REAL, PARAMETER :: KARMAN = 0.4
	150	!...
	151
	152	RDT = 1.0 / DTPBL
	153
	154	DO K = 1, kte
	155	SIGMAF(K-1) = ZNW(K)
	156	ENDDO
	157	SIGMAF(kte) = 0.0
	158
	159	DO K = 1, kte
	160	DSIGH(K) = SIGMAF(K) - SIGMAF(K-1)
	161	DSIGHI(K) = 1.0 / DSIGH(K)
	162	ENDDO
	163
	164	DO K = kts,kte-1
	165	DSIGFI(K) = 1.0 / (SIGMAH(K+1) - SIGMAH(K))
	166	ENDDO
	167
	168	DSIGFI(kte) = DSIGFI(kte-1)
	169	do j = jts,jte
	170	CALL ACM2D(j=J,xtime=XTIME, dtpbl=DTPBL, sigmaf=SIGMAF, sigmah=SIGMAH &
	171	,dsigfi=DSIGFI,dsighi=DSIGHI,dsigh=DSIGH &
	172	,us=u3d(ims,kms,j),vs=v3d(ims,kms,j) &
	173	,theta=th3d(ims,kms,j),tt=t3d(ims,kms,j) &
	174	,qvs=qv3d(ims,kms,j),qcs=qc3d(ims,kms,j) &
	175	,qis=qi3d(ims,kms,j),dzf=DZ8W(ims,kms,j) &
	176	,densx=RR3D(ims,kms,j) &
	177	,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) &
	178	,ttnp=rthblten(ims,kms,j),qvtnp=rqvblten(ims,kms,j) &
	179	,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &
	180	,cpd=cpd,g=g,rovcp=rovcp,rd=rd,rdt=rdt &
	181	,psfcpa=psfc(ims,j),ust=ust(ims,j) &
	182	,pbl=pblh(ims,j) &
	183	,regime=regime(ims,j),psim=psim(ims,j) &
	184	,hfx=hfx(ims,j),qfx=qfx(ims,j) &
	185	,tg=tsk(ims,j),gz1oz0=gz1oz0(ims,j) &
	186	,wspd=wspd(ims,j) ,klpbl=kpbl2d(ims,j) &
	187	,mut=mut(ims,j) &
	188	,ep1=ep1,karman=karman &
	189	,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
	190	,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
	191	,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
	192	enddo
	193
	194	END SUBROUTINE ACMPBL
	195
	196	!--------------------------------------------------------------------
	197	SUBROUTINE ACM2D(j,XTIME, DTPBL, sigmaf, sigmah &
	198	,dsigfi,dsighi,dsigh &
	199	,us,vs,theta,tt,qvs,qcs,qis &
	200	,dzf,densx,utnp,vtnp,ttnp,qvtnp,qctnp,qitnp &
	201	,cpd,g,rovcp,rd,rdt,psfcpa,ust &
	202	,pbl,regime,psim &
	203	,hfx,qfx,tg,gz1oz0,wspd ,klpbl &
	204	,mut &
	205	,ep1,karman &
	206	,ids,ide, jds,jde, kds,kde &
	207	,ims,ime, jms,jme, kms,kme &
	208	,its,ite, jts,jte, kts,kte )
	209
	210	IMPLICIT NONE
	211
	212	!.......Arguments
	213
	214	!... Real
	215	REAL, DIMENSION( 0:kte ), INTENT(IN) :: SIGMAF
	216	REAL, DIMENSION( kms:kme ), INTENT(IN) :: SIGMAH
	217	REAL, DIMENSION( kts:kte ), INTENT(IN) :: DSIGH, DSIGHI, DSIGFI
	218	REAL , INTENT(IN) :: DTPBL, G, RD,ep1,karman,CPD,ROVCP,RDT
	219	REAL , DIMENSION( ims:ime ), INTENT(INOUT) :: PBL, UST
	220
	221	REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, TT, &
	222	QVS, QCS, QIS, DENSX
	223	REAL, DIMENSION( ims:ime, kms:kme ), intent(in) :: DZF
	224	REAL, DIMENSION( ims:ime, kms:kme ), intent(inout) :: utnp, &
	225	vtnp, &
	226	ttnp, &
	227	qvtnp, &
	228	qctnp, &
	229	qitnp
	230	real, dimension( ims:ime ), intent(in ) :: psfcpa
	231	real, dimension( ims:ime ), intent(in ) :: tg
	232	real, dimension( ims:ime ), intent(inout) :: regime
	233	real, dimension( ims:ime ), intent(in) :: wspd, psim, gz1oz0
	234	real, dimension( ims:ime ), intent(in) :: hfx, qfx
	235	real, dimension( ims:ime ), intent(in) :: mut
	236
	237	!... Integer
	238	INTEGER, DIMENSION( ims:ime ), INTENT(OUT):: KLPBL
	239	INTEGER, INTENT(IN) :: XTIME
	240	integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
	241	ims,ime, jms,jme, kms,kme, &
	242	its,ite, jts,jte, kts,kte, j
	243	!--------------------------------------------------------------------
	244	!--Local
	245	INTEGER I, K
	246	INTEGER :: KPBLHT
	247	INTEGER, DIMENSION( its:ite ) :: KPBLH, NOCONV
	248
	249	!... Real
	250	REAL :: TVCON, WSS, TCONV, TH1, TOG, DTMP, WSSQ
	251	REAL :: psix, THV1
	252	REAL, DIMENSION( its:ite ) :: FINT, PSTAR, CPAIR
	253	REAL, DIMENSION( its:ite, kts:kte ) :: THETAV, RIB, &
	254	EDDYZ, UX, VX, THETAX, &
	255	QVX, QCX, QIX, ZA
	256	REAL, DIMENSION( its:ite, 0:kte ) :: ZF
	257	REAL, DIMENSION( its:ite) :: WST, TST, QST, USTM, TSTV
	258	REAL, DIMENSION( its:ite ) :: PBLSIG, MOL
	259	REAL :: FINTT, ZMIX, UMIX, VMIX
	260	REAL :: TMPFX, TMPVTCON, TMPP, TMPTHS, TMPTH1, TMPVCONV, WS1, DTH
	261	REAL :: A,TST12,RL,ZFUNC
	262	! REAL, PARAMETER :: KARMAN = 0.4
	263
	264	!... Integer
	265	INTEGER :: KL, jtf, ktf, itf, KMIX
	266	!...
	267	!-----initialize vertical tendencies and
	268	!
	269	do i = its,ite
	270	do k = kts,kte
	271	utnp(i,k) = 0.
	272	vtnp(i,k) = 0.
	273	ttnp(i,k) = 0.
	274	enddo
	275	enddo
	276	!
	277	do k = kts,kte
	278	do i = its,ite
	279	qvtnp(i,k) = 0.
	280	enddo
	281	enddo
	282	!
	283	do k = kts,kte
	284	do i = its,ite
	285	qctnp(i,k) = 0.
	286	qitnp(i,k) = 0.
	287	enddo
	288	enddo
	289	!
	290	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	291	!!! Compute Micromet Scaling variables, not availiable in WRF for ACM
	292	DO I = its,ite
	293	cpair(i) = CPD * (1.0 + 0.84 * QVS(I,1)) ! J/(K KG)
	294	TMPFX = HFX(I) / (cpair(i) * DENSX(I,1))
	295	TMPVTCON = 1.0 + EP1 * QVS(I,1) ! COnversion factor for virtual temperature
	296	WS1 = SQRT(US(I,1)2 + VS(I,1)2) ! Level 1 wind speed
	297	TST(I) = -TMPFX / UST(I)
	298	QST(I) = -QFX(I) / (UST(I)*DENSX(I,1))
	299	USTM(I) = UST(I) * WS1 / wspd(i)
	300	THV1 = TMPVTCON * THETA(I,1)
	301	TSTV(I) = TST(I)TMPVTCON + THV1EP1*QST(I)
	302	MOL(I) = THV1 * UST(i)*2/(KARMANG*TSTV(I))
	303	WST(I) = UST(I) * (PBL(I)/(KARMANABS(MOL(I)))) * 0.333333
	304	PSTAR(I) = MUT(I)/1000. ! P* in cb
	305	ENDDO
	306	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	307
	308	!... Compute PBL height
	309
	310	!... compute the height of full- and half-sigma level above ground level
	311	DO I = its,ite
	312	ZF(I,0) = 0.0
	313	KLPBL(I) = 1
	314	ENDDO
	315
	316	DO K = kts, kte
	317	DO I = its,ite
	318	ZF(I,K) = DZF(I,K) + ZF(I,K-1)
	319	ZA(I,K) = 0.5 * (ZF(I,K) + ZF(I,K-1))
	320	ENDDO
	321	ENDDO
	322
	323	DO K = kts, kte
	324	DO I = its,ite
	325	TVCON = 1.0 + EP1 * QVS(I,K)
	326	THETAV(I,K) = THETA(I,K) * TVCON
	327	ENDDO
	328	ENDDO
	329
	330
	331	!... COMPUTE PBL WHERE RICHARDSON NUMBER = RIC (0.25) HOLTSLAG ET AL 1990
	332	DO 100 I = its,ite
	333	if(MOL(I).LT.0.0) then
	334	WSS = (UST(I) ** 3 + 0.6 * WST(I) 3) 0.33333
	335	TMPFX = -UST(I) * TSTV(I)
	336	TCONV = 8.5 * TMPFX / WSS
	337	TH1 = THETAV(I,1) + TCONV
	338	else
	339	TH1 = THETAV(I,1)
	340	endif
	341	99 KMIX = 1
	342	DO K = 1,kte
	343	DTMP = THETAV(I,K) - TH1
	344	IF (DTMP.LT.0.0) KMIX = K
	345	ENDDO
	346	IF(KMIX.GT.1) THEN
	347	FINTT = (TH1 - THETAV(I,KMIX)) / (THETAV(I,KMIX+1) &
	348	- THETAV(I,KMIX))
	349	ZMIX = FINTT * (ZA(I,KMIX+1)-ZA(I,KMIX)) + ZA(I,KMIX)
	350	UMIX = FINTT * (US(I,KMIX+1)-US(I,KMIX)) + US(I,KMIX)
	351	VMIX = FINTT * (VS(I,KMIX+1)-VS(I,KMIX)) + VS(I,KMIX)
	352	ELSE
	353	ZMIX = ZA(I,1)
	354	UMIX = US(I,1)
	355	VMIX = VS(I,1)
	356	ENDIF
	357	DO K = KMIX,kte
	358	DTMP = THETAV(I,K) - TH1
	359	TOG = 0.5 * (THETAV(I,K) + TH1) / G
	360	WSSQ = (US(I,K)-UMIX)**2 &
	361	+ (VS(I,K)-VMIX)**2
	362	IF (KMIX == 1) WSSQ = WSSQ + 100.UST(I)UST(I)
	363	WSSQ = MAX( WSSQ, 0.1 )
	364	RIB(I,K) = ABS(ZA(I,K)-ZMIX) * DTMP / (TOG * WSSQ)
	365	IF (RIB(I,K) .GE. RIC) GO TO 201
	366	ENDDO
	367
	368	print *,' RIB never exceeds RIC, RIB(i,kte) = ',rib(i,5), &
	369	' THETAV(i,1) = ',thetav(i,1),' MOL=',mol(i), &
	370	' TCONV = ',TCONV,' WST = ',WST(I), &
	371	' KMIX = ',kmix,' UST = ',UST(I), &
	372	' TST = ',TST(I),' U,V = ',US(I,1),VS(I,1), &
	373	' I,J=',I,J
	374	STOP
	375	201 CONTINUE
	376
	377	KPBLH(I) = K
	378
	379	100 CONTINUE
	380
	381	DO I = its,ite
	382	IF (KPBLH(I) .NE. 1) THEN
	383	!---------INTERPOLATE BETWEEN LEVELS -- jp 7/93
	384	FINT(I) = (RIC - RIB(I,KPBLH(I)-1)) / (RIB(I,KPBLH(I)) - &
	385	RIB(I,KPBLH(I)-1))
	386	IF (FINT(I) .GT. 0.5) THEN
	387	KPBLHT = KPBLH(I)
	388	FINT(I) = FINT(I) - 0.5
	389	ELSE
	390	KPBLHT = KPBLH(I) - 1
	391	FINT(I) = FINT(I) + 0.5
	392	ENDIF
	393	PBL(I) = FINT(I) * (ZF(I,KPBLHT) - ZF(I,KPBLHT-1)) + &
	394	ZF(I,KPBLHT-1)
	395	KLPBL(I) = KPBLHT
	396	PBLSIG(I) = FINT(I) * DSIGH(KPBLHT) + SIGMAF(KPBLHT-1) ! sigma at PBL height
	397	ELSE
	398	KLPBL(I) = 1
	399	PBL(I) = ZF(I,1)
	400	PBLSIG(I) = SIGMAF(1)
	401	ENDIF
	402
	403	ENDDO
	404
	405	DO I = its,ite
	406	NOCONV(I) = 0
	407
	408	! Check for CBL and identify conv. vs. non-conv cells
	409	IF (PBL(I) / MOL(I) .LT. -0.02 .AND. KLPBL(I) .GT. 3 &
	410	.AND. THETAV(I,1) .GT. THETAV(I,2) .AND. XTIME .GT. 1) THEN
	411	NOCONV(I) = 1
	412	REGIME(I) = 4.0 ! FREE CONVECTIVE - ACM
	413	ENDIF
	414	ENDDO
	415
	416	!... Calculate Kz
	417	CALL EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, &
	418	US, VS, TT, THETAV, DENSX, PSTAR, &
	419	QVS, QCS, QIS, DSIGFI, G, RD, CPAIR, &
	420	EDDYZ, its,ite, kts,kte,ims,ime, kms,kme)
	421
	422	CALL ACM(DTPBL, PSTAR, NOCONV, SIGMAF, DSIGH, DSIGHI, J, &
	423	KLPBL, PBL, PBLSIG, MOL, UST, &
	424	TST, QST, USTM, EDDYZ, DENSX, &
	425	US, VS, THETA, QVS, QCS, QIS, &
	426	UX, VX, THETAX, QVX, QCX, QIX, &
	427	ids,ide, jds,jde, kds,kde, &
	428	ims,ime, jms,jme, kms,kme, &
	429	its,ite, jts,jte, kts,kte)
	430
	431	!... Calculate tendency due to PBL parameterization
	432
	433	DO K = kts, kte
	434	DO I = its, ite
	435	UTNP(I,K) = UTNP(I,K) + (UX(I,K) - US(I,K)) * RDT
	436	VTNP(I,K) = VTNP(I,K) + (VX(I,K) - VS(I,K)) * RDT
	437	TTNP(I,K) = TTNP(I,K) + (THETAX(I,K) - THETA(I,K)) * RDT
	438	QVTNP(I,K) = QVTNP(I,K) + (QVX(I,K) - QVS(I,K)) * RDT
	439	QCTNP(I,K) = QCTNP(I,K) + (QCX(I,K) - QCS(I,K)) * RDT
	440	QITNP(I,K) = QITNP(I,K) + (QIX(I,K) - QIS(I,K)) * RDT
	441	ENDDO
	442	ENDDO
	443	! IF(J.EQ.36)PRINT *,' PBL,THETA,THETAX,QVS,QVX=',PBL(20),THETA(20,1),THETAX(20,1),QVS(20,1),QVX(20,1)
	444	! IF(J.EQ.36)PRINT *,' UST,ustm,TG,TST=',UST(20),ustm(20),Tg(20),TST(20)
	445	! IF(J.EQ.36)PRINT *,' HFX,MOL,WST,WSPD=',HFX(20),MOL(20),WST(20),wspd(20)
	446	! IF(J.EQ.36)then
	447	! i=20
	448	! do k=1,kte
	449	! PRINT *,' qvten,uten,vten=',QVTNP(I,K),UTNP(I,K),VTNP(I,K)
	450	! PRINT *,' k,thten,th,thx,edy=',k,TTNP(I,K),THETA(20,k),THETAX(20,k),EDDYZ(20,K)
	451	! enddo
	452	! ENDIF
	453	END SUBROUTINE ACM2D
	454
	455	!-------------------------------------------------------------------
	456	SUBROUTINE ACMINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
	457	RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &
	458	restart, allowed_to_read , &
	459	ids, ide, jds, jde, kds, kde, &
	460	ims, ime, jms, jme, kms, kme, &
	461	its, ite, jts, jte, kts, kte )
	462	!**********************************************************************
	463	!
	464	! This subroutine is for preparing ACM PBL variables.
	465	! Called from module_physics_init.F
	466	!
	467	! REVISION HISTORY:
	468	! AX 3/2005 - Originally developed
	469	!**********************************************************************
	470	! ARGUMENT LIST:
	471	!
	472	!-------------------------------------------------------------------
	473	IMPLICIT NONE
	474	!
	475	LOGICAL , INTENT(IN) :: restart , allowed_to_read
	476
	477	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
	478	ims, ime, jms, jme, kms, kme, &
	479	its, ite, jts, jte, kts, kte
	480
	481	INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR
	482
	483	! REAL , DIMENSION( kms:kme ), INTENT(IN) :: SHALF
	484	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
	485	RUBLTEN, &
	486	RVBLTEN, &
	487	RTHBLTEN, &
	488	RQVBLTEN, &
	489	RQCBLTEN, &
	490	RQIBLTEN
	491
	492	!... Local Variables
	493	INTEGER :: i, j, k, itf, jtf, ktf
	494
	495	!
	496	jtf=min0(jte,jde-1)
	497	ktf=min0(kte,kde-1)
	498	itf=min0(ite,ide-1)
	499
	500	IF(.not.restart)THEN
	501	DO j=jts,jtf
	502	DO k=kts,ktf
	503	DO i=its,itf
	504	RUBLTEN(i,k,j)=0.
	505	RVBLTEN(i,k,j)=0.
	506	RTHBLTEN(i,k,j)=0.
	507	RQVBLTEN(i,k,j)=0.
	508	RQCBLTEN(i,k,j)=0.
	509	ENDDO
	510	ENDDO
	511	ENDDO
	512	ENDIF
	513
	514	IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
	515	DO j=jts,jtf
	516	DO k=kts,ktf
	517	DO i=its,itf
	518	RQIBLTEN(i,k,j)=0.
	519	ENDDO
	520	ENDDO
	521	ENDDO
	522	ENDIF
	523
	524
	525	END SUBROUTINE acminit
	526
	527	!-------------------------------------------------------------------
	528	SUBROUTINE EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, &
	529	US, VS, TT, THETAV, DENSX, PSTAR, &
	530	QVS, QCS, QIS, DSIGFI, G, RD, CPAIR, &
	531	EDDYZ, its,ite, kts,kte,ims,ime,kms,kme )
	532
	533
	534	!**********************************************************************
	535	! Two methods for computing Kz:
	536	! 1. Boundary scaling similar to Holtslag and Boville (1993)
	537	! 2. Local Kz computed as function of local Richardson # and vertical
	538	! wind shear, similar to LIU & CARROLL (1996)
	539	!
	540	!**********************************************************************
	541	!
	542	!-- DTPBL time step of the minor loop for the land-surface/pbl model
	543	!-- ZF height of full sigma level
	544	!-- ZA height of half sigma level
	545	!-- MOL Monin-Obukhov length in 1D form
	546	!-- PBL PBL height in 1D form
	547	!-- UST friction velocity U* in 1D form (m/s)
	548	!-- US U wind
	549	!-- VS V wind
	550	!-- TT temperature
	551	!-- THETAV potential virtual temperature
	552	!-- DENSX dry air density (kg/m^3)
	553	!-- PSTAR P*=Psfc-Ptop
	554	!-- QVS water vapor mixing ratio (Kg/Kg)
	555	!-- QCS cloud mixing ratio (Kg/Kg)
	556	!-- QIS ice mixing ratio (Kg/Kg)
	557	!-- DSIGFI inverse of sigma layer delta
	558	!-- G gravity
	559	!-- RD gas constant for dry air (j/kg/k)
	560	!-- CPAIR specific heat of moist air (M^2 S^-2 K^-1)
	561	!-- EDDYZ eddy diffusivity KZ
	562	!-----------------------------------------------------------------------
	563
	564	IMPLICIT NONE
	565
	566	!.......Arguments
	567
	568	!... Integer
	569	INTEGER, INTENT(IN) :: its,ite, kts,kte,ims,ime,kms,kme
	570	!... Real
	571	REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST
	572	REAL , INTENT(IN) :: DTPBL, G, RD
	573	REAL , DIMENSION( kts:kte ), INTENT(IN) :: DSIGFI
	574	REAL , DIMENSION( its:ite ), INTENT(IN) :: MOL, PSTAR, CPAIR
	575
	576	REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, TT, &
	577	QVS, QCS, QIS, DENSX
	578	REAL, DIMENSION( its:ite, kts:kte ), INTENT(IN) :: ZA, THETAV
	579	REAL, DIMENSION( its:ite, 0:kte ) , INTENT(IN) :: ZF
	580
	581	REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: EDDYZ
	582
	583	!.......Local variables
	584
	585	!... Integer
	586	INTEGER :: ILX, KL, KLM, K, I
	587
	588	!... Real
	589	REAL :: ZOVL, PHIH, WT, ZSOL, ZFUNC, DZF, SS, GOTH, EDYZ
	590	REAL :: RI, QMEAN, TMEAN, XLV, ALPH, CHI, ZK, SQL, DENSF,kzo
	591
	592	!... Parameters
	593	REAL, PARAMETER :: RV = 461.5
	594	REAL, PARAMETER :: RC = 0.25
	595	REAL, PARAMETER :: RLAM = 80.0
	596	REAL, PARAMETER :: GAMH = 16.0 !15.0 ! Holtslag and Boville (1993)
	597	REAL, PARAMETER :: BETAH = 5.0 ! Holtslag and Boville (1993)
	598	REAL, PARAMETER :: KARMAN = 0.4
	599	REAL, PARAMETER :: EDYZ0 = 0.01 ! New Min Kz
	600	! REAL, PARAMETER :: EDYZ0 = 0.1
	601	!-- IMVDIF imvdif=1 for moist adiabat vertical diffusion
	602	INTEGER, PARAMETER :: imvdif = 1
	603	!
	604	ILX = ite
	605	KL = kte
	606	KLM = kte - 1
	607
	608	DO K = kts,KLM
	609	DO I = its,ILX
	610	EDYZ = 0.0
	611	ZOVL = 0.0
	612	DZF = ZA(I,K+1) - ZA(I,K)
	613	! kzo = min(0.001 * dzf,EDYZ0)
	614	kzo = 0.001 * dzf
	615	! kzo = EDYZ0
	616	!-------------------------------------------------
	617	IF (ZF(I,K) .LT. PBL(I)) THEN
	618	ZOVL = ZF(I,K) / MOL(I)
	619	IF (ZOVL .LT. 0.0) THEN
	620	IF (ZF(I,K) .LT. 0.1 * PBL(I)) THEN
	621	PHIH = 1.0 / SQRT(1.0 - GAMH * ZOVL)
	622	WT = UST(I) / PHIH
	623	ELSE
	624	ZSOL = 0.1 * PBL(I) / MOL(I)
	625	PHIH = 1.0 / SQRT(1.0 - GAMH * ZSOL)
	626	WT = UST(I) / PHIH
	627	ENDIF
	628	ELSE IF (ZOVL .LT. 1.0) THEN
	629	PHIH = 1.0 + BETAH * ZOVL
	630	WT = UST(I) / PHIH
	631	ELSE
	632	PHIH = BETAH + ZOVL
	633	WT = UST(I) / PHIH
	634	ENDIF
	635	ZFUNC = ZF(I,K) * (1.0 - ZF(I,K) / PBL(I)) ** 2
	636	EDYZ = KARMAN * WT * ZFUNC
	637	! EDYZ = AMAX1(EDYZ,EDYZ0)
	638	ENDIF
	639	! PRINT *,I,K,TT(I,K)
	640	!--------------------------------------------------------------------------
	641	!
	642	! RC = 0.257 * DZF ** 0.175
	643	SS = ((US(I,K+1) - US(I,K)) 2 + (VS(I,K+1) - VS(I,K)) 2) &
	644	/ (DZF * DZF) + 1.0E-9
	645	GOTH = 2.0 * G / (THETAV(I,K+1) + THETAV(I,K))
	646	RI = GOTH * (THETAV(I,K+1) - THETAV(I,K)) / (DZF * SS)
	647	! PRINT *,I,K,TT(I,K), RI, IMVDIF
	648	!
	649	!-- Adjustment to vert diff in Moist air
	650	if(imvdif.eq.1)then
	651	IF ((QCS(I,K)+QIS(I,K)) .GT. 0.01E-3 .OR. (QCS(I,K+1)+ &
	652	QIS(I,K+1)) .GT. 0.01E-3) THEN
	653	QMEAN = 0.5 * (QVS(I,K) + QVS(I,K+1))
	654	TMEAN = 0.5 * (TT(I,K) + TT(I,K+1))
	655	XLV = (2.501 - 0.00237 * (TMEAN - 273.15)) * 1.E6
	656	ALPH = XLV * QMEAN / RD / TMEAN
	657	CHI = XLV * XLV * QMEAN / CPAIR(I) / RV / TMEAN / TMEAN
	658	RI = (1.0 + ALPH) * (RI -G * G / SS / TMEAN / CPAIR(I) * &
	659	((CHI - ALPH) / (1.0 + CHI)))
	660	ENDIF
	661	endif
	662
	663	ZK = 0.4 * ZF(I,K)
	664	SQL = (ZK * RLAM / (RLAM + ZK)) ** 2
	665
	666	IF (RI .GE. RC) THEN
	667	EDDYZ(I,K) = kzo !EDYZ0
	668	ELSE IF (RI .GE. 0.0) THEN
	669	EDDYZ(I,K) = kzo + SQRT(SS) * (1.- RI / RC) ** 2 * SQL
	670	ELSE
	671	EDDYZ(I,K) = kzo + SQRT(SS * (1.0 - 25.0 * RI)) * SQL
	672	ENDIF
	673
	674	IF(EDYZ.GT.EDDYZ(I,K)) THEN
	675	EDDYZ(I,K) = EDYZ
	676	ENDIF
	677
	678	EDDYZ(I,K) = MIN(1000.0,EDDYZ(I,K))
	679	EDDYZ(I,K) = MAX(kzo,EDDYZ(I,K))
	680
	681	DENSF = 0.5 * (DENSX(I,K+1) + DENSX(I,K))
	682
	683	EDDYZ(I,K) = EDDYZ(I,K) * (DENSF * G / PSTAR(I)) ** 2 * &
	684	DTPBL * DSIGFI(K)*1E-6
	685	ENDDO ! for I loop
	686	ENDDO ! for k loop
	687	!
	688	DO I = its,ILX
	689	EDDYZ(I,KL) = 0.0 ! EDDYZ(I,KLM) -- changed jp 3/08
	690	ENDDO
	691	END SUBROUTINE EDDYX
	692
	693	!-------------------------------------------------------------------
	694	SUBROUTINE ACM (DTPBL, PSTAR, NOCONV, SIGMAF, DSIGH, DSIGHI, JX, &
	695	KLPBL, PBL, PBLSIG, MOL, UST, &
	696	TST, QST, USTM, EDDYZ, DENSX, &
	697	US, VS, THETA, QVS, QCS, QIS, &
	698	UX, VX, THETAX, QVX, QCX, QIX, &
	699	ids,ide, jds,jde, kds,kde, &
	700	ims,ime, jms,jme, kms,kme, &
	701	its,ite, jts,jte, kts,kte)
	702	!**********************************************************************
	703	! PBL model called the Asymmetric Convective Model, Version 2 (ACM2)
	704	! -- See top of module for summary and references
	705	!
	706	!---- REVISION HISTORY:
	707	! AX 3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM
	708	! JP and RG 8/2006 - updates
	709	!
	710	!**********************************************************************
	711	! ARGUMENTS:
	712	!-- DTPBL PBL time step
	713	!-- PSTAR Psurf - Ptop in cb
	714	!-- NOCONV If free convection =0, no; =1, yes
	715	!-- SIGMAF Sigma for full layer
	716	!-- DSIGH Sigma thickness
	717	!-- DSIGHI Inverse of sigma thickness
	718	!-- JX N-S index
	719	!-- KLPBL PBL level at K index
	720	!-- PBL PBL height in m
	721	!-- PBLSIG Sigma level for PBL
	722	!-- MOL Monin-Obukhov length in 1D form
	723	!-- UST U* in 1D form
	724	!-- TST Theta* in 1D form
	725	!-- QST Q* in 1D form
	726	!-- USTM U* for computation of momemtum flux
	727	!-- EDDYZ eddy diffusivity KZ
	728	!-- DENSX dry air density (kg/m^3)
	729	!-- US U wind
	730	!-- VS V wind
	731	!-- THETA potential temperature
	732	!-- QVS water vapor mixing ratio (Kg/Kg)
	733	!-- QCS cloud mixing ratio (Kg/Kg)
	734	!-- QIS ice mixing ratio (Kg/Kg)
	735	!-- UX new U wind
	736	!-- VX new V wind
	737	!-- THETAX new potential temperature
	738	!-- QVX new water vapor mixing ratio (Kg/Kg)
	739	!-- QCX new cloud mixing ratio (Kg/Kg)
	740	!-- QIX new ice mixing ratio (Kg/Kg)
	741	!-----------------------------------------------------------------------
	742
	743	IMPLICIT NONE
	744
	745	!.......Arguments
	746
	747	!... Integer
	748	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
	749	ims,ime, jms,jme, kms,kme, &
	750	its,ite, jts,jte, kts,kte, JX
	751	INTEGER, DIMENSION( its:ite ), INTENT(IN) :: NOCONV
	752	INTEGER, DIMENSION( ims:ime ), INTENT(IN) :: KLPBL
	753
	754	!... Real
	755	REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST
	756	REAL , INTENT(IN) :: DTPBL
	757	REAL , DIMENSION( its:ite ), INTENT(IN) :: PSTAR, PBLSIG, &
	758	MOL, TST, &
	759	QST, USTM
	760	REAL , DIMENSION( kts:kte ), INTENT(IN) :: DSIGHI, DSIGH
	761	REAL , DIMENSION( 0:kte ), INTENT(IN) :: SIGMAF
	762	REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT) :: EDDYZ
	763	REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, &
	764	QVS, QCS, QIS, DENSX
	765	REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: UX, VX, THETAX, &
	766	QVX, QCX, QIX
	767
	768	!.......Local variables
	769
	770	!... Parameters
	771	INTEGER, PARAMETER :: NSP = 6
	772	!
	773	!......ACM2 Parameters
	774	! INTEGER, PARAMETER :: IFACM = 0
	775	!
	776	REAL, PARAMETER :: G1000 = 9.8 * 1.0E-3
	777	REAL, PARAMETER :: XX = 0.5 ! FACTOR APPLIED TO CONV MIXING TIME STEP
	778	REAL, PARAMETER :: KARMAN = 0.4
	779
	780	!... Integer
	781	INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L
	782	INTEGER :: KCBLMX
	783	INTEGER, DIMENSION( its:ite ) :: KCBL
	784
	785	!... Real
	786	REAL :: G1000I, MBMAX, HOVL, MEDDY, MBAR
	787	REAL :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1
	788	REAL, DIMENSION( its:ite ) :: PSTARI, FSACM, RAH, DTLIM
	789	REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN
	790	REAL, DIMENSION( 1:NSP, its:ite ) :: FS, BCBOTN
	791	REAL, DIMENSION( kts:kte ) :: XPLUS, XMINUS
	792	REAL DELC
	793	REAL, DIMENSION( 1:NSP,its:ite,kts:kte ) :: VCI
	794
	795	REAL, DIMENSION( kts:kte ) :: AI, BI, CI, EI !, Y
	796	REAL, DIMENSION( 1:NSP,kts:kte ) :: DI, UI
	797	!
	798	!--Start Exicutable ----
	799
	800	ILX = ite
	801	KL = kte
	802	KLM = kte - 1
	803
	804	G1000I = 1.0 / G1000
	805	KCBLMX = 0
	806	MBMAX = 0.0
	807
	808	!---COMPUTE ACM MIXING RATE
	809	DO I = its, ILX
	810	DTLIM(I) = DTPBL
	811	PSTARI(I) = 1.0 / PSTAR(I)
	812	KCBL(I) = 1
	813	FSACM(I) = 0.0
	814
	815	IF (NOCONV(I) .EQ. 1) THEN
	816	KCBL(I) = KLPBL(I)
	817	!-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE
	818	!--New couple ACM & EDDY-------------------------------------------------------------
	819	HOVL = -PBL(I) / MOL(I)
	820	FSACM(I) = 1./(1.+((KARMAN/(HOVL))*0.3333)/(0.72KARMAN))
	821	MEDDY = EDDYZ(I,1) / (DTPBL * (PBLSIG(I) - SIGMAF(1)))
	822	MBAR = MEDDY * FSACM(I)
	823	DO K = kts,KCBL(I)-1
	824	EDDYZ(I,K) = EDDYZ(I,K) * (1.0 - FSACM(I))
	825	ENDDO
	826	! if(i.eq.100.and.jx.eq.43) PRINT *,' Edy,meddy,mbar=', EDDYZ(I,1),MEDDY,MBAR
	827	MBMAX = AMAX1(MBMAX,MBAR)
	828	DO K = kts+1,KCBL(I)
	829	MBARKS(K,I) = MBAR
	830	MDWN(K,I) = MBAR * (PBLSIG(I) - SIGMAF(K-1)) * DSIGHI(K)
	831	ENDDO
	832	MBARKS(1,I) = MBAR
	833	MBARKS(KCBL(I),I) = MDWN(KCBL(I),I)
	834	MDWN(KCBL(I)+1,I) = 0.0
	835	ENDIF
	836	ENDDO ! end of I loop
	837
	838	! if(NOCONV(20).EQ.1) print *,' KCBL,PBLSIG=',KCBL(20),PBLSIG(20),SIGMAF(KCBL(20)-1),DSIGHI(KCBL(20))
	839	! do k = kts,klm
	840	! if(NOCONV(20).EQ.1.and.MBAR.lt.1.0e-3)print *,' k,MBARKS,MDWN=',k,MBARKS(k,20),MDWN(k,20)
	841	! enddo
	842
	843	DO K = kts,KLM
	844	DO I = its,ILX
	845	EKZ = EDDYZ(I,K) / DTPBL * DSIGHI(K)
	846	DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I))
	847	ENDDO
	848	ENDDO
	849
	850	DO I = its,ILX
	851	IF (NOCONV(I) .EQ. 1) THEN
	852	KCBLMX = AMAX0(KLPBL(I),KCBLMX)
	853	RZ = (SIGMAF(KCBL(I)) - SIGMAF(1)) * DSIGHI(1)
	854	DTLIM(I) = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I))
	855	ENDIF
	856	ENDDO
	857
	858	DO I = its,ILX
	859	ENDDO
	860	!
	861	DO K = kts,KL
	862	DO I = its,ILX
	863	VCI(1,I,K) = THETA(I,K)
	864	VCI(2,I,K) = QVS(I,K)
	865	VCI(3,I,K) = US(I,K)
	866	VCI(4,I,K) = VS(I,K)
	867	! -- Also mix cloud water and ice if necessary
	868	! IF (IMOISTX.NE.1.AND.IMOISTX.NE.3) THEN !!! Check other PBL models
	869	VCI(5,I,K) = QCS(I,K)
	870	VCI(6,I,K) = QIS(I,K)
	871	ENDDO
	872	ENDDO
	873
	874	NSPX=6
	875
	876	DO I = its,ILX
	877	! RAH(I) = RA(I,J) + 3.976 / UST(I)
	878	FS(1,I) = -UST(I) * TST(I) * DENSX(I,1) * PSTARI(I)
	879	FS(2,I) = -UST(I) * QST(I) * DENSX(I,1) * PSTARI(I)
	880	FM = -USTM(I) * USTM(I) * DENSX(I,1) * PSTARI(I)
	881	WSPD = SQRT(US(I,1) * US(I,1) + VS(I,1) * VS(I,1)) + 1.E-9
	882	FS(3,I) = FM * US(I,1) / WSPD
	883	FS(4,I) = FM * VS(I,1) / WSPD
	884	FS(5,I) = 0.0
	885	FS(6,I) = 0.0 ! SURFACE FLUXES OF CLOUD WATER AND ICE = 0
	886	ENDDO
	887	!
	888	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	889	DO I = its,ILX
	890	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	891	NLP = INT(DTPBL / DTLIM(I) + 1.0)
	892	DTS = (DTPBL / NLP)
	893	DTRAT = DTS / DTPBL
	894	DO NL = 1,NLP ! LOOP OVER SUB TIME LOOP
	895	!
	896	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	897	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	898
	899	!-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES
	900
	901	DO K = kts,KL
	902	AI(K) = 0.0
	903	BI(K) = 0.0
	904	CI(K) = 0.0
	905	EI(K) = 0.0
	906	ENDDO
	907
	908	DO K = 2, KCBL(I)
	909	EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DSIGH(K) * DSIGHI(K-1)
	910	BI(K) = 1.0 + CRANKP * MDWN(K,I) * DTS
	911	AI(K) = -CRANKP * MBARKS(K,I) * DTS
	912	ENDDO
	913
	914	EI(1) = EI(1) -EDDYZ(I,1) * CRANKP * DSIGHI(1 )* DTRAT
	915	AI(2) = AI(2) -EDDYZ(I,1) * CRANKP * DSIGHI(2) * DTRAT
	916
	917	DO K = KCBL(I)+1, KL
	918	BI(K) = 1.0
	919	ENDDO
	920
	921	DO K = 2,KL
	922	XPLUS(K) = EDDYZ(I,K) * DSIGHI(K) * DTRAT
	923	XMINUS(K) = EDDYZ(I,K-1) * DSIGHI(K) * DTRAT
	924	CI(K) = - XMINUS(K) * CRANKP
	925	EI(K) = EI(K) - XPLUS(K) * CRANKP
	926	BI(K) = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP
	927	ENDDO
	928
	929	IF (NOCONV(I) .EQ. 1) THEN
	930	BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBLSIG(I) - SIGMAF(1)) * &
	931	DTS * DSIGHI(1) + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
	932	ELSE
	933	BI(1) = 1.0 + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
	934	ENDIF
	935
	936
	937	DO K = 1,KL
	938	DO L = 1,NSPX
	939	DI(L,K) = 0.0
	940	ENDDO
	941	ENDDO
	942	!
	943	!** COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION
	944	DO K = 2,KCBL(I)
	945	DO L = 1,NSPX
	946	DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) * &
	947	VCI(L,I,K) + DSIGH(K+1) * DSIGHI(K) * &
	948	MDWN(K+1,I) * VCI(L,I,K+1))
	949	DI(L,K) = VCI(L,I,K) + (1.0 - CRANKP) * DELC
	950	ENDDO
	951	ENDDO
	952
	953	DO K = KCBL(I)+1, KL
	954	DO L = 1,NSPX
	955	DI(L,K) = VCI(L,I,K)
	956	ENDDO
	957	ENDDO
	958
	959	DO K = 2,KL
	960	IF (K .EQ. KL) THEN
	961	DO L = 1,NSPX
	962	DI(L,K) = DI(L,K) - (1.0 - CRANKP) * XMINUS(K) * &
	963	(VCI(L,I,K) - VCI(L,I,K-1))
	964	ENDDO
	965	ELSE
	966	DO L = 1,NSPX
	967	DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) * &
	968	(VCI(L,I,K+1) - VCI(L,I,K)) - &
	969	(1.0 - CRANKP) * XMINUS(K) * &
	970	(VCI(L,I,K) - VCI(L,I,K-1))
	971	ENDDO
	972	ENDIF
	973	ENDDO
	974
	975	IF (NOCONV(I) .EQ. 1) THEN
	976	DO L = 1,NSPX
	977	F1 = -G1000I * (MBARKS(1,I) * &
	978	(PBLSIG(I) - SIGMAF(1)) * VCI(L,I,1) - &
	979	MDWN(2,I) * VCI(L,I,2) * DSIGH(2))
	980
	981	DI(L,1) = VCI(L,I,1) - G1000 * (FS(L,I) - (1.0 - CRANKP) &
	982	* F1) * DSIGHI(1) * DTS
	983	ENDDO
	984	ELSE
	985	DO L = 1,NSPX
	986	DI(L,1) = VCI(L,I,1) - G1000 * FS(L,I) * DSIGHI(1) * DTS
	987	ENDDO
	988	ENDIF
	989	DO L = 1,NSPX
	990	DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZ(I,1) * DSIGHI(1) &
	991	* DTRAT * (VCI(L,I,2) - VCI(L,I,1))
	992	ENDDO
	993	IF ( NOCONV(I) .EQ. 1 ) THEN
	994	CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX)
	995	ELSE
	996	CALL TRI (CI, BI, EI, DI, UI, KL, NSPX)
	997	END IF
	998	!
	999	!-- COMPUTE NEW THETAV AND Q
	1000	DO K = 1,KL
	1001	DO L = 1,NSPX
	1002	VCI(L,I,K) = UI(L,K)
	1003	ENDDO
	1004	ENDDO
	1005	!
	1006	ENDDO ! END I LOOP
	1007	ENDDO ! END SUB TIME LOOP
	1008	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1009	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1010
	1011	!
	1012	DO K = kts,KL
	1013	DO I = its,ILX
	1014	THETAX(I,K) = VCI(1,I,K)
	1015	QVX(I,K) = VCI(2,I,K)
	1016	UX(I,K) = VCI(3,I,K)
	1017	VX(I,K) = VCI(4,I,K)
	1018	! IF( (I.EQ.61 .OR. I.EQ.15) .AND. JX.EQ.22) THEN
	1019	! PRINT *,'TEST -->',I,K,THETAX(I,K)
	1020	! ENDIF
	1021	ENDDO
	1022	ENDDO
	1023
	1024	DO K = kts,KL
	1025	DO I = its,ILX
	1026	QCX(I,K) = VCI(5,I,K)
	1027	QIX(I,K) = VCI(6,I,K)
	1028	ENDDO
	1029	ENDDO
	1030	! ENDIF
	1031
	1032	END SUBROUTINE ACM
	1033	!--------------------------------------------------------
	1034
	1035	SUBROUTINE MATRIX(A,B,C,D,E,X,KL,NSP)
	1036
	1037	IMPLICIT NONE
	1038	!
	1039	!-- Bordered band diagonal matrix solver for ACM2
	1040
	1041	!-- ACM2 Matrix is in this form:
	1042	! B1 E1
	1043	! A2 B2 E2
	1044	! A3 C3 B3 E3
	1045	! A4 C4 B4 E4
	1046	! A5 C5 B5 E5
	1047	! A6 C6 B6
	1048
	1049	!--Upper Matrix is
	1050	! U11 U12
	1051	! U22 U23
	1052	! U33 U34
	1053	! U44 U45
	1054	! U55 U56
	1055	! U66
	1056
	1057	!--Lower Matrix is:
	1058	! 1
	1059	! L21 1
	1060	! L31 L32 1
	1061	! L41 L42 L43 1
	1062	! L51 L52 L53 L54 1
	1063	! L61 L62 L63 L64 L65 1
	1064	!---------------------------------------------------------
	1065	!...Arguments
	1066	INTEGER, INTENT(IN) :: KL
	1067	INTEGER, INTENT(IN) :: NSP
	1068	REAL A(KL),B(KL),E(KL)
	1069	REAL C(KL),D(NSP,KL),X(NSP,KL)
	1070
	1071	!...Locals
	1072	REAL Y(NSP,KL),AIJ,SUM
	1073	REAL L(KL,KL),UII(KL),UIIP1(KL),RUII(KL)
	1074	INTEGER I,J,V
	1075
	1076	!-- Define Upper and Lower matrices
	1077	L(1,1) = 1.
	1078	! U(1,1) = B(1)
	1079	UII(1) = B(1)
	1080	RUII(1) = 1./UII(1)
	1081	DO I = 2, KL
	1082	L(I,I) = 1.
	1083	L(I,1) = A(I)/B(1)
	1084	! U(I-1,I) =E(I-1)
	1085	UIIP1(I-1)=E(I-1)
	1086	IF(I.GE.3) THEN
	1087	DO J = 2,I-1
	1088	IF(I.EQ.J+1) THEN
	1089	AIJ = C(I)
	1090	ELSE
	1091	AIJ = 0.
	1092	ENDIF
	1093	L(I,J) = (AIJ-L(I,J-1)*E(J-1))/ &
	1094	(B(J)-L(J,J-1)*E(J-1))
	1095	ENDDO
	1096	ENDIF
	1097	ENDDO
	1098
	1099	DO I = 2,KL
	1100	! U(I,I) = B(I)-L(I,I-1)*E(I-1)
	1101	UII(I) = B(I)-L(I,I-1)*E(I-1)
	1102	RUII(I) = 1./UII(I)
	1103	ENDDO
	1104
	1105	!-- Forward sub for Ly=d
	1106	DO V= 1, NSP
	1107	Y(V,1) = D(V,1)
	1108	DO I=2,KL
	1109	SUM = D(V,I)
	1110	DO J=1,I-1
	1111	SUM = SUM - L(I,J)*Y(V,J)
	1112	ENDDO
	1113	Y(V,I) = SUM
	1114	ENDDO
	1115	ENDDO
	1116
	1117	!-- Back sub for Ux=y
	1118
	1119	DO V= 1, NSP
	1120	! X(V,KL) = Y(V,KL)/U(KL,KL)
	1121	X(V,KL) = Y(V,KL)*RUII(KL)
	1122	ENDDO
	1123	DO I = KL-1,1,-1
	1124	DO V= 1, NSP
	1125	! X(V,I) = (Y(V,I)-U(I,I+1)*X(V,I+1))/U(I,I)
	1126	X(V,I) = (Y(V,I)-UIIP1(I)X(V,I+1))RUII(I)
	1127	ENDDO
	1128	ENDDO
	1129
	1130	END SUBROUTINE MATRIX
	1131	!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
	1132	SUBROUTINE TRI ( L, D, U, B, X,KL,NSP)
	1133	!-----------------------------------------------------------------------
	1134
	1135	! FUNCTION:
	1136	! Solves tridiagonal system by Thomas algorithm.
	1137	! The associated tri-diagonal system is stored in 3 arrays
	1138	! D : diagonal
	1139	! L : sub-diagonal
	1140	! U : super-diagonal
	1141	! B : right hand side function
	1142	! X : return solution from tridiagonal solver
	1143
	1144	! [ D(1) U(1) 0 0 0 ... 0 ]
	1145	! [ L(2) D(2) U(2) 0 0 ... . ]
	1146	! [ 0 L(3) D(3) U(3) 0 ... . ]
	1147	! [ . . . . . ] X(i) = B(i)
	1148	! [ . . . . 0 ]
	1149	! [ . . . . ]
	1150	! [ 0 L(n) D(n) ]
	1151
	1152	!-----------------------------------------------------------------------
	1153
	1154	IMPLICIT NONE
	1155
	1156	! Arguments:
	1157
	1158	INTEGER, INTENT(IN) :: KL
	1159	INTEGER, INTENT(IN) :: NSP
	1160
	1161	REAL L( KL ) ! subdiagonal
	1162	REAL D(KL) ! diagonal
	1163	REAL U( KL ) ! superdiagonal
	1164	REAL B(NSP,KL ) ! R.H. side
	1165	REAL X( NSP,KL ) ! solution
	1166
	1167	! Local Variables:
	1168
	1169	REAL GAM( KL )
	1170	REAL BET
	1171	INTEGER V, K
	1172
	1173	! Decomposition and forward substitution:
	1174	BET = 1.0 / D( 1 )
	1175	DO V = 1, NSP
	1176	X( V,1 ) = BET * B(V,1 )
	1177	ENDDO
	1178
	1179	DO K = 2, KL
	1180	GAM(K ) = BET * U( K-1 )
	1181	BET = 1.0 / ( D( K ) - L( K ) * GAM( K ) )
	1182	DO V = 1, NSP
	1183	X( V, K ) = BET * ( B( V,K ) - L( K ) * X( V,K-1 ) )
	1184	ENDDO
	1185	END DO
	1186
	1187	! Back-substitution:
	1188
	1189	DO K = KL - 1, 1, -1
	1190	DO V = 1, NSP
	1191	X( V,K ) = X( V,K ) - GAM( K+1 ) * X( V,K+1 )
	1192	END DO
	1193	ENDDO
	1194
	1195	END SUBROUTINE TRI
	1196
	1197	END MODULE module_bl_acm
	1198

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

Download in other formats:

Original Format