Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_bl_mrf.F @ 3553

Last change on this file since 3553 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 80.1 KB

Rev	Line
[2759]	1	!WRF:MODEL_LAYER:PHYSICS
	2	!
	3	MODULE module_bl_mrf
	4
	5	CONTAINS
	6
	7	!-------------------------------------------------------------------
	8	SUBROUTINE MRF(U3D,V3D,TH3D,T3D,QV3D,QC3D,P3D,PI3D, &
	9	RUBLTEN,RVBLTEN,RTHBLTEN, &
	10	RQVBLTEN,RQCBLTEN, &
	11	CP,G,ROVCP,R,ROVG, &
	12	dz8w,z,XLV,RV,PSFC, &
	13	p1000mb, &
	14	ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH, &
	15	XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR, &
	16	DT,DTMIN,KPBL2D, &
	17	SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
	18	flag_qi, &
	19	ids,ide, jds,jde, kds,kde, &
	20	ims,ime, jms,jme, kms,kme, &
	21	its,ite, jts,jte, kts,kte, &
	22	! Optional
	23	QI3D,RQIBLTEN, &
	24	regime )
	25	!-------------------------------------------------------------------
	26	IMPLICIT NONE
	27	!-------------------------------------------------------------------
	28	!-- U3D 3D u-velocity interpolated to theta points (m/s)
	29	!-- V3D 3D v-velocity interpolated to theta points (m/s)
	30	!-- TH3D 3D potential temperature (K)
	31	!-- T3D temperature (K)
	32	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
	33	!-- QC3D 3D cloud mixing ratio (Kg/Kg)
	34	!-- QI3D 3D ice mixing ratio (Kg/Kg)
	35	!-- P3D 3D pressure (Pa)
	36	!-- PI3D 3D exner function (dimensionless)
	37	!-- rr3D 3D dry air density (kg/m^3)
	38	!-- RUBLTEN U tendency due to
	39	! PBL parameterization (m/s^2)
	40	!-- RVBLTEN V tendency due to
	41	! PBL parameterization (m/s^2)
	42	!-- RTHBLTEN Theta tendency due to
	43	! PBL parameterization (K/s)
	44	!-- RQVBLTEN Qv tendency due to
	45	! PBL parameterization (kg/kg/s)
	46	!-- RQCBLTEN Qc tendency due to
	47	! PBL parameterization (kg/kg/s)
	48	!-- RQIBLTEN Qi tendency due to
	49	! PBL parameterization (kg/kg/s)
	50	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
	51	!-- G acceleration due to gravity (m/s^2)
	52	!-- ROVCP R/CP
	53	!-- R gas constant for dry air (J/kg/K)
	54	!-- ROVG R/G
	55	!-- dz8w dz between full levels (m)
	56	!-- z height above sea level (m)
	57	!-- XLV latent heat of vaporization (J/kg)
	58	!-- RV gas constant for water vapor (J/kg/K)
	59	!-- PSFC pressure at the surface (Pa)
	60	!-- ZNT roughness length (m)
	61	!-- UST u* in similarity theory (m/s)
	62	!-- ZOL z/L height over Monin-Obukhov length
	63	!-- HOL PBL height over Monin-Obukhov length
	64	!-- PBL PBL height (m)
	65	!-- REGIME flag indicating PBL regime (stable, unstable, etc.)
	66	!-- PSIM similarity stability function for momentum
	67	!-- PSIH similarity stability function for heat
	68	!-- XLAND land mask (1 for land, 2 for water)
	69	!-- HFX upward heat flux at the surface (W/m^2)
	70	!-- QFX upward moisture flux at the surface (kg/m^2/s)
	71	!-- TSK surface temperature (K)
	72	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
	73	!-- WSPD wind speed at lowest model level (m/s)
	74	!-- BR bulk Richardson number in surface layer
	75	!-- DT time step (s)
	76	!-- DTMIN time step (minute)
	77	!-- rvovrd R_v divided by R_d (dimensionless)
	78	!-- SVP1 constant for saturation vapor pressure (kPa)
	79	!-- SVP2 constant for saturation vapor pressure (dimensionless)
	80	!-- SVP3 constant for saturation vapor pressure (K)
	81	!-- SVPT0 constant for saturation vapor pressure (K)
	82	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
	83	!-- EP2 constant for specific humidity calculation
	84	!-- KARMAN Von Karman constant
	85	!-- EOMEG angular velocity of earth's rotation (rad/s)
	86	!-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
	87	!-- ids start index for i in domain
	88	!-- ide end index for i in domain
	89	!-- jds start index for j in domain
	90	!-- jde end index for j in domain
	91	!-- kds start index for k in domain
	92	!-- kde end index for k in domain
	93	!-- ims start index for i in memory
	94	!-- ime end index for i in memory
	95	!-- jms start index for j in memory
	96	!-- jme end index for j in memory
	97	!-- kms start index for k in memory
	98	!-- kme end index for k in memory
	99	!-- its start index for i in tile
	100	!-- ite end index for i in tile
	101	!-- jts start index for j in tile
	102	!-- jte end index for j in tile
	103	!-- kts start index for k in tile
	104	!-- kte end index for k in tile
	105	!-------------------------------------------------------------------
	106
	107	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	108	ims,ime, jms,jme, kms,kme, &
	109	its,ite, jts,jte, kts,kte
	110
	111	!
	112	REAL, INTENT(IN ) :: P1000mb
	113	REAL, INTENT(IN ) :: DT,DTMIN,CP,G,ROVCP, &
	114	ROVG,R,XLV,RV
	115
	116	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
	117	REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
	118
	119	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	120	INTENT(IN ) :: QV3D, &
	121	QC3D, &
	122	P3D, &
	123	PI3D, &
	124	TH3D, &
	125	T3D, &
	126	dz8w, &
	127	z
	128	!
	129	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	130	INTENT(INOUT) :: RUBLTEN, &
	131	RVBLTEN, &
	132	RTHBLTEN, &
	133	RQVBLTEN, &
	134	RQCBLTEN
	135
	136	REAL, DIMENSION( ims:ime, jms:jme ) , &
	137	INTENT(IN ) :: XLAND, &
	138	HFX, &
	139	QFX
	140
	141	REAL, DIMENSION( ims:ime, jms:jme ) , &
	142	INTENT(INOUT) :: HOL, &
	143	UST, &
	144	PBL, &
	145	ZNT
	146
	147	LOGICAL, INTENT(IN) :: FLAG_QI
	148	!
	149	!m The following 5 variables are changed to memory size from tile size--
	150	!
	151	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: &
	152	PSIM, &
	153	PSIH
	154
	155	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
	156	WSPD
	157
	158	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: &
	159	GZ1OZ0, &
	160	BR
	161
	162	REAL, DIMENSION( ims:ime, jms:jme ) , &
	163	INTENT(IN ) :: PSFC
	164
	165	REAL, DIMENSION( ims:ime, jms:jme ) , &
	166	INTENT(IN ) :: TSK
	167
	168	REAL, DIMENSION( ims:ime, jms:jme ) , &
	169	INTENT(INOUT) :: ZOL
	170
	171	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
	172	INTENT(OUT ) :: KPBL2D
	173
	174	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	175	INTENT(IN ) :: U3D, &
	176	V3D
	177	!
	178	! Optional
	179	!
	180	REAL, DIMENSION( ims:ime, jms:jme ) , &
	181	OPTIONAL , &
	182	INTENT(INOUT) :: REGIME
	183
	184	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	185	INTENT(INOUT) :: RQIBLTEN
	186
	187	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	188	OPTIONAL , &
	189	INTENT(IN ) :: QI3D
	190
	191	! LOCAL VARS
	192	REAL, DIMENSION( its:ite, kts:kte ) :: dz8w2d, &
	193	z2d
	194
	195
	196	INTEGER :: I,J,K,NK
	197
	198	!
	199	DO J=jts,jte
	200	DO k=kts,kte
	201	NK=kme-k
	202	DO i=its,ite
	203	dz8w2d(I,K) = dz8w(i,NK,j)
	204	z2d(I,K) = z(i,NK,j)
	205	ENDDO
	206	ENDDO
	207
	208
	209	CALL MRF2D(J,U3D(ims,kms,j),V3D(ims,kms,j),T3D(ims,kms,j), &
	210	QV3D(ims,kms,j),QC3D(ims,kms,j), &
	211	P3D(ims,kms,j),RUBLTEN(ims,kms,j),RVBLTEN(ims,kms,j),&
	212	RTHBLTEN(ims,kms,j),RQVBLTEN(ims,kms,j), &
	213	RQCBLTEN(ims,kms,j), &
	214	p1000mb, &
	215	CP,G,ROVCP,R,ROVG, &
	216	dz8w2d,z2d,XLV,Rv, &
	217	PSFC(ims,j),ZNT(ims,j), &
	218	UST(ims,j),ZOL(ims,j), &
	219	HOL(ims,j),PBL(ims,j),PSIM(ims,j), &
	220	PSIH(ims,j),XLAND(ims,j),HFX(ims,j),QFX(ims,j), &
	221	TSK(ims,j),GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j), &
	222	DT,DTMIN,KPBL2D(ims,j), &
	223	SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT, &
	224	flag_qi, &
	225	ids,ide, jds,jde, kds,kde, &
	226	ims,ime, jms,jme, kms,kme, &
	227	its,ite, jts,jte, kts,kte, &
	228	!optional
	229	QI2DTEN=RQIBLTEN(ims,kms,j), &
	230	REGIME=REGIME(ims,j),QI2D=QI3D(ims,kms,j) )
	231
	232
	233	DO k=kts,kte
	234	DO i=its,ite
	235	RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/PI3D(I,K,J)
	236	ENDDO
	237	ENDDO
	238
	239	ENDDO
	240
	241	END SUBROUTINE MRF
	242
	243	!-------------------------------------------------------------------
	244	SUBROUTINE MRF2D(J,U2D,V2D,T2D,QV2D,QC2D, P2D, &
	245	U2DTEN,V2DTEN,T2DTEN, &
	246	QV2DTEN,QC2DTEN, &
	247	p1000mb, &
	248	CP,G,ROVCP,R,ROVG, &
	249	dz8w2d,z2d,XLV,RV,PSFCPA, &
	250	ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH, &
	251	XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR, &
	252	DT,DTMIN,KPBL1D, &
	253	SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
	254	flag_qi, &
	255	ids,ide, jds,jde, kds,kde, &
	256	ims,ime, jms,jme, kms,kme, &
	257	its,ite, jts,jte, kts,kte, &
	258	! optional
	259	regime, qi2d, QI2DTEN )
	260	!-------------------------------------------------------------------
	261	IMPLICIT NONE
	262	!-------------------------------------------------------------------
	263	! BASED ON THE "COUNTERGRADIENT" TRANSPORT TERM OF TROEN
	264	! AND MAHRT (1986) FOR THE UNSTABLE PBL.
	265	! THIS ROUTINE USES AN IMPLICIT APPROACH FOR VERTICAL FLUX
	266	! DIVERGENCE AND DOES NOT REQUIRE "MITER" TIMESTEPS.
	267	! IT INCLUDES VERTICAL DIFFUSION IN THE STABLE ATMOSPHERE
	268	! AND MOIST VERTICAL DIFFUSION IN CLOUDS.
	269	! SURFACE FLUXES CALCULATED AS IN HIRPBL.
	270	! 5-LAYER SOIL MODEL OPTION REQUIRED IN SLAB DUE TO LONG TIMESTEP
	271	!
	272	! CODED BY SONG-YOU HONG (NCEP), IMPLEMENTED BY JIMY DUDHIA (NCAR)
	273	! FALL 1996
	274	!
	275	! REFERENCES:
	276	!
	277	! HONG AND PAN (1996), MON. WEA. REV.
	278	! TROEN AND MAHRT (1986), BOUNDARY LAYER MET.
	279	!
	280	! CHANGES:
	281	! INCREASE RLAM FROM 30 TO 150, AND CHANGE FREE ATMOSPHERE
	282	! STABILITY FUNCTION TO INCREASE VERTICAL DIFFUSION
	283	! (HONG, JUNE 1997)
	284	!
	285	! PUT LOWER LIMIT ON PSI FOR STABLE CONDITIONS. THIS WILL
	286	! PREVENT FLUXES FROM BECOMING TOO SMALL (DUDHIA, OCTOBER 1997)
	287	!
	288	! CORRECTION TO REGIME CALCULATION. THIS WILL ALLOW POINTS IN
	289	! REGIME 4 MUCH MORE FREQUENTLY GIVING LARGER SURFACE FLUXES
	290	! REGIME 3 NO LONGER USES HOL < 1.5 OR THVX LAPSE-RATE CHECK
	291	! IN MRF SCHEME. THIS WILL MAKE REGIME 3 MUCH LESS FREQUENT.
	292	!
	293	! ADD SURFACE PRESSURE, PS(I), ARRAY FOR EFFICIENCY
	294	!
	295	! FIX FOR PROBLEM WITH THIN LAYERS AND HIGH ROUGHNESS
	296	!
	297	! CHARNOCK CONSTANT NOW COMES FROM NAMELIST (DEFAULT SAME)
	298	!
	299	!-------------------------------------------------------------------
	300
	301	REAL RLAM,PRMIN,PRMAX,XKZMIN,XKZMAX,RIMIN,BRCR, &
	302	CFAC,PFAC,SFCFRAC,CKZ,ZFMIN,APHI5,APHI16,GAMCRT, &
	303	GAMCRQ,XKA,PRT
	304
	305	PARAMETER (RLAM=150.,PRMIN=0.5,PRMAX=4.)
	306	PARAMETER (XKZMIN=0.01,XKZMAX=1000.,RIMIN=-100.)
	307	PARAMETER (BRCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1)
	308	PARAMETER (CKZ=0.001,ZFMIN=1.E-8,APHI5=5.,APHI16=16.)
	309	PARAMETER (GAMCRT=3.,GAMCRQ=2.E-3)
	310	PARAMETER (XKA=2.4E-5)
	311	PARAMETER (PRT=1.)
	312	!
	313	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	314	ims,ime, jms,jme, kms,kme, &
	315	its,ite, jts,jte, kts,kte, &
	316	J
	317	!
	318	LOGICAL, INTENT(IN) :: FLAG_QI
	319	!
	320	REAL, INTENT(IN ) :: P1000mb
	321	REAL, INTENT(IN ) :: DT,DTMIN,CP,G,ROVCP, &
	322	ROVG,R,XLV,RV
	323
	324	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
	325	REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
	326
	327	REAL, DIMENSION( ims:ime, kms:kme ) , &
	328	INTENT(IN ) :: QV2D, &
	329	QC2D, &
	330	P2D, &
	331	T2D
	332	!
	333	REAL, DIMENSION( ims:ime, kms:kme ) , &
	334	INTENT(INOUT) :: U2DTEN, &
	335	V2DTEN, &
	336	T2DTEN, &
	337	QV2DTEN, &
	338	QC2DTEN
	339
	340	REAL, DIMENSION( ims:ime ) , &
	341	INTENT(INOUT) :: HOL, &
	342	UST, &
	343	PBL, &
	344	ZNT
	345
	346	REAL, DIMENSION( ims:ime ) , &
	347	INTENT(IN ) :: XLAND, &
	348	HFX, &
	349	QFX
	350	!
	351	!m The following 5 are changed to memory size---
	352	!
	353	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: PSIM, &
	354	PSIH
	355
	356	REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: WSPD
	357
	358	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: GZ1OZ0, &
	359	BR
	360
	361	REAL, DIMENSION( ims:ime ) , &
	362	INTENT(IN ) :: PSFCPA
	363
	364	REAL, DIMENSION( ims:ime ) , &
	365	INTENT(IN ) :: TSK
	366
	367	REAL, DIMENSION( ims:ime ) , &
	368	INTENT(INOUT) :: ZOL
	369
	370	INTEGER, DIMENSION( ims:ime ) , &
	371	INTENT(OUT ) :: KPBL1D
	372
	373	REAL, DIMENSION( ims:ime, kms:kme ) , &
	374	INTENT(IN ) :: U2D, &
	375	V2D
	376
	377	! MODULE-LOCAL VARIABLES (DEFINED IN SUBROUTINE MRF)
	378	!
	379	REAL, DIMENSION( its:ite, kts:kte ) , &
	380	INTENT(IN) :: dz8w2d, &
	381	z2d
	382	!
	383	!
	384	! Optional
	385	!
	386	REAL, DIMENSION( ims:ime ) , &
	387	OPTIONAL , &
	388	INTENT(INOUT) :: REGIME
	389
	390	REAL, DIMENSION( ims:ime, kms:kme ) , &
	391	OPTIONAL , &
	392	INTENT(IN ) :: QI2D
	393
	394	REAL, DIMENSION( ims:ime, kms:kme ) , &
	395	OPTIONAL , &
	396	INTENT(INOUT) :: QI2DTEN
	397
	398	! LOCAL VARS
	399
	400	REAL, DIMENSION( its:ite, kts:kte+1 ) :: ZQ
	401
	402	REAL, DIMENSION( its:ite, kts:kte ) :: &
	403	UX,VX,QX, &
	404	QCX,THX,THVX, &
	405	DZQ,DZA, &
	406	TTNP,QTNP, &
	407	QCTNP,ZA, &
	408	UXS,VXS, &
	409	THXS,QXS, &
	410	QCXS,QIX, &
	411	QITNP,QIXS, &
	412	UTNP,VTNP
	413	!
	414	REAL, DIMENSION( its:ite ) :: QIXSV,RHOX, &
	415	WSPD1,GOVRTH, &
	416	PBL0,THXSV, &
	417	UXSV,VXSV, &
	418	QXSV,QCXSV, &
	419	QGH,TGDSA,PS
	420
	421	INTEGER :: ILXM,JLXM,KL, &
	422	KLM,KLP1,KLPBL
	423	!
	424	INTEGER, DIMENSION( its:ite ) :: KPBL,KPBL0
	425	!
	426	REAL, DIMENSION( its:ite, kts:kte ) :: SCR3,SCR4
	427	!
	428	REAL, DIMENSION( its:ite ) :: DUM1, &
	429	XKZMKL
	430	!
	431	REAL, DIMENSION( its:ite ) :: ZL1,THERMAL, &
	432	WSCALE,HGAMT, &
	433	HGAMQ,BRDN, &
	434	BRUP,PHIM, &
	435	PHIH,CPM, &
	436	DUSFC,DVSFC, &
	437	DTSFC,DQSFC
	438
	439	!
	440	REAL, DIMENSION( its:ite, kts:kte ) :: XKZM,XKZH, &
	441	A1,A2, &
	442	AD,AU, &
	443	TX
	444	!
	445	REAL, DIMENSION( its:ite, kts:kte ) :: AL
	446	!
	447	LOGICAL, DIMENSION( its:ite ) :: PBLFLG, &
	448	SFCFLG, &
	449	STABLE
	450	!
	451	REAL, DIMENSION( its:ite ) :: THGB
	452
	453	INTEGER :: N,I,K,KK,L,NZOL,IMVDIF
	454
	455	INTEGER :: JBGN,JEND,IBGN,IEND,NK
	456
	457	REAL :: ZOLN,X,Y,CONT,CONQ,CONW,PL,THCON,TVCON,E1,DTSTEP
	458	REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL
	459	REAL :: DTTHX,PSIX,DTG,PSIQ,USTM
	460	REAL :: DT4,RDT,SPDK2,FM,FH,HOL1,GAMFAC,VPERT,PRNUM
	461	REAL :: ZFAC,XKZO,SS,RI,QMEAN,TMEAN,ALPH,CHI,ZK,RL2,DK,SRI
	462	REAL :: BRINT,DTODSD,DSIG,RDZ,DSDZT,DSDZQ,DSDZ2,TTEND,QTEND
	463	REAL :: UTEND,VTEND,QCTEND,QITEND,TGC,DTODSU
	464
	465	!----------------------------------------------------------------------
	466
	467	KLPBL=1
	468	KL=kte
	469	ILXM=ite-1
	470	JLXM=jte-1
	471	KLM=kte-1
	472	KLP1=kte+1
	473	!
	474	CONT=1000.*CP/G
	475	CONQ=1000.*XLV/G
	476	CONW=1000./G
	477
	478	!-- IMVDIF imvdif=1 for moist adiabat vertical diffusion
	479
	480	IMVDIF=1
	481
	482	! DO i=its,ite
	483	!!PS PSFC cmb
	484	! PSFC(I)=PSFCPA(I)/1000.
	485	! ENDDO
	486
	487
	488	!
	489	!----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
	490	!
	491	DO 5 I=its,ite
	492	TGDSA(I)=TSK(I)
	493	! PS PSFC cmb
	494	PS(I)=PSFCPA(I)/1000.
	495	! THGB(I)=TSK(I)(100./PS(I))*ROVCP
	496	THGB(I)=TSK(I)(P1000mb/PSFCPA(I))*ROVCP
	497	5 CONTINUE
	498	!
	499	!-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
	500	! T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
	501	!
	502	! * NOTE *
	503	! THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
	504	! SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
	505	! TENDENCIES.
	506	!
	507	DO 24 K=kts,kte
	508	NK=kme-K
	509	DO 24 I=its,ite
	510	UX(I,K)=U2D(I,NK)
	511	VX(I,K)=V2D(I,NK)
	512	24 CONTINUE
	513	!
	514	!.....SCR3(I,K) STORE TEMPERATURE,
	515	! SCR4(I,K) STORE VIRTUAL TEMPERATURE.
	516	!
	517	DO 30 K=kts,kte
	518	NK=kme-K
	519	DO 30 I=its,ite
	520	! PL cmb
	521	PL=P2D(I,NK)/1000.
	522	SCR3(I,K)=T2D(I,NK)
	523	! THCON=(100./PL)**ROVCP
	524	THCON=(P1000mb/(PL1000.))*ROVCP
	525	THX(I,K)=SCR3(I,K)*THCON
	526	TX(I,K)=SCR3(I,K)
	527	SCR4(I,K)=SCR3(I,K)
	528	THVX(I,K)=THX(I,K)
	529	QX(I,K)=0.
	530	30 CONTINUE
	531	!
	532	DO I=its,ite
	533	QGH(i)=0.
	534	CPM(i)=CP
	535	ENDDO
	536	!
	537	! IF(IDRY.EQ.1)GOTO 80
	538	DO 50 K=kts,kte
	539	NK=kme-K
	540	DO 50 I=its,ite
	541	QX(I,K)=QV2D(I,NK)
	542	TVCON=(1.+EP1*QX(I,K))
	543	THVX(I,K)=THX(I,K)*TVCON
	544	SCR4(I,K)=SCR3(I,K)*TVCON
	545	50 CONTINUE
	546	!
	547	DO 60 I=its,ite
	548	E1=SVP1EXP(SVP2(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
	549	QGH(I)=EP2*E1/(PS(I)-E1)
	550	CPM(I)=CP(1.+0.8QX(I,KL))
	551	60 CONTINUE
	552	!
	553	! IF(IMOIST.EQ.1)GOTO 80
	554	DO 70 K=kts,kte
	555	NK=kme-K
	556	DO 70 I=its,ite
	557	QCX(I,K)=QC2D(I,NK)
	558	70 CONTINUE
	559
	560	IF (flag_QI .AND. PRESENT( QI2D ) ) THEN
	561	DO K=kts,kte
	562	NK=kme-K
	563	DO I=its,ite
	564	QIX(I,K)=QI2D(I,NK)
	565	ENDDO
	566	ENDDO
	567	ELSE
	568	DO K=kts,kte
	569	NK=kme-K
	570	DO I=its,ite
	571	QIX(I,K)=0.
	572	ENDDO
	573	ENDDO
	574	ENDIF
	575
	576	80 CONTINUE
	577
	578	!
	579	!-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
	580	! LEVEL, AND THE LAYER THICKNESSES.
	581	!
	582	DO 90 I=its,ite
	583	ZQ(I,KLP1)=0.
	584	RHOX(I)=PS(I)1000./(RSCR4(I,KL))
	585	90 CONTINUE
	586	!
	587	DO 110 KK=kts,kte
	588	K=kme-KK
	589	DO 100 I=its,ite
	590	DUM1(I)=ZQ(I,K+1)
	591	100 CONTINUE
	592	!
	593	DO 110 I=its,ite
	594	ZQ(I,K)=dz8w2d(I,K)+DUM1(I)
	595	110 CONTINUE
	596	!
	597	DO 120 K=kts,kte
	598	DO 120 I=its,ite
	599	ZA(I,K)=0.5*(ZQ(I,K)+ZQ(I,K+1))
	600	DZQ(I,K)=ZQ(I,K)-ZQ(I,K+1)
	601	120 CONTINUE
	602	!
	603	DO 130 K=kts,kte-1
	604	DO 130 I=its,ite
	605	DZA(I,K)=ZA(I,K)-ZA(I,K+1)
	606	130 CONTINUE
	607
	608	DTSTEP=DT
	609	!
	610	DO 160 I=its,ite
	611	GOVRTH(I)=G/THX(I,KL)
	612	160 CONTINUE
	613	!
	614	!-----INITIALIZE VERTICAL TENDENCIES AND
	615	!
	616	DO I=its,ite
	617	DO K=kts,kte
	618	UTNP(i,k)=0.
	619	VTNP(i,k)=0.
	620	TTNP(i,k)=0.
	621	ENDDO
	622	ENDDO
	623	!
	624	! IF(IDRY.EQ.1)GOTO 250
	625	DO 230 K=kts,kte
	626	DO 230 I=its,ite
	627	QTNP(I,K)=0.
	628	230 CONTINUE
	629	!
	630	! IF(IMOIST.EQ.1)GOTO 250
	631	DO 240 K=kts,kte
	632	DO 240 I=its,ite
	633	QCTNP(I,K)=0.
	634	QITNP(I,K)=0.
	635	240 CONTINUE
	636
	637	250 CONTINUE
	638	!
	639	!-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
	640	! AKB(1976), EQ(12).
	641
	642	! DO 260 I=its,ite
	643	! GZ1OZ0(I)=ALOG(ZA(I,KL)/ZNT(I))
	644	! IF((XLAND(I)-1.5).GE.0)THEN
	645	! ZL=ZNT(I)
	646	! ELSE
	647	! ZL=0.01
	648	! ENDIF
	649	! WSPD(I)=SQRT(UX(I,KL)UX(I,KL)+VX(I,KL)VX(I,KL))
	650	! TSKV=THGB(I)(1.+EP1QGH(I)*MAVAIL(I))
	651	! DTHVDZ=(THVX(I,KL)-TSKV)
	652	! IF(-DTHVDZ.GE.0)THEN
	653	! DTHVM=-DTHVDZ
	654	! ELSE
	655	! DTHVM=0.
	656	! ENDIF
	657	! VCONV=VCONVC*SQRT(DTHVM)
	658	! WSPD(I)=SQRT(WSPD(I)WSPD(I)+VCONVVCONV)
	659	! WSPD(I)=AMAX1(WSPD(I),1.)
	660	! BR(I)=GOVRTH(I)ZA(I,KL)DTHVDZ/(WSPD(I)*WSPD(I))
	661	! 260 CONTINUE
	662
	663	!!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
	664	!!
	665	!! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
	666	!! AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
	667	!!
	668	!! CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
	669	!!
	670	!! 1. BR .GE. 0.2;
	671	!! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
	672	!!
	673	!! 2. BR .LT. 0.2 .AND. BR .GT. 0.0;
	674	!! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
	675	!! (REGIME=2),
	676	!!
	677	!! 3. BR .EQ. 0.0
	678	!! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
	679	!!
	680	!! 4. BR .LT. 0.0
	681	!! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
	682	!!
	683	!!-----
	684	!
	685	! DO 320 I=its,ite
	686	!!----
	687	!!-- REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
	688	!!-- IF(BR(I).LT.0..AND.HOL(I).GT.1.5)GOTO 310
	689	!
	690	! IF(BR(I).LT.0.)GOTO 310
	691	!!
	692	!!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
	693	!!
	694	! IF(BR(I).LT.0.2)GOTO 270
	695	! REGIME(I)=1.
	696	! PSIM(I)=-10.*GZ1OZ0(I)
	697	!! LOWER LIMIT ON PSI IN STABLE CONDITIONS
	698	! PSIM(I)=AMAX1(PSIM(I),-10.)
	699	! PSIH(I)=PSIM(I)
	700	! HOL(I)=0.0
	701	! PBL(I)=0.0
	702	! GOTO 320
	703	!!
	704	!!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
	705	!!
	706	! 270 IF(BR(I).EQ.0.0)GOTO 280
	707	! REGIME(I)=2.
	708	! PSIM(I)=-5.0BR(I)GZ1OZ0(I)/(1.1-5.0*BR(I))
	709	!! LOWER LIMIT ON PSI IN STABLE CONDITIONS
	710	! PSIM(I)=AMAX1(PSIM(I),-10.)
	711	!!.....AKB(1976), EQ(16).
	712	! PSIH(I)=PSIM(I)
	713	! HOL(I)=0.0
	714	! PBL(I)=0.0
	715	! GOTO 320
	716	!!
	717	!!-----CLASS 3; FORCED CONVECTION:
	718	!!
	719	! 280 REGIME(I)=3.
	720	! PSIM(I)=0.0
	721	! PSIH(I)=PSIM(I)
	722	!
	723	!! special use kte instead of kme
	724	!
	725	! DO 290 KK=kts,kte-1
	726	! K=kte-KK
	727	! IF(THVX(I,K).GT.THVX(I,KL))GOTO 300
	728	! 290 CONTINUE
	729	! STOP 290
	730	! 300 PBL(I)=ZQ(I,K+1)
	731	! IF(UST(I).LT.0.01)THEN
	732	! ZOL(I)=BR(I)*GZ1OZ0(I)
	733	! ELSE
	734	! ZOL(I)=KARMANGOVRTH(I)ZA(I,KL)MOL(I,J)/(UST(I)UST(I))
	735	! ENDIF
	736	! HOL(I)=-ZOL(I)*PBL(I)/ZA(I,KL)
	737	! GOTO 320
	738	!
	739	!!-----CLASS 4; FREE CONVECTION:
	740	!
	741	!! 310 IF(THVX(I,KLM).GT.THVX(I,KL))GOTO 280
	742	!
	743	! 310 CONTINUE
	744	! REGIME(I)=4.
	745	! IF(UST(I).LT.0.01)THEN
	746	! ZOL(I)=BR(I)*GZ1OZ0(I)
	747	! ELSE
	748	! ZOL(I)=KARMANGOVRTH(I)ZA(I,KL)MOL(I,J)/(UST(I)UST(I))
	749	! ENDIF
	750	! ZOL(I)=AMIN1(ZOL(I),0.)
	751	! ZOL(I)=AMAX1(ZOL(I),-9.9999)
	752	! NZOL=INT(-ZOL(I)*100.)
	753	! RZOL=-ZOL(I)*100.-NZOL
	754	! PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
	755	! PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
	756	!!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
	757	!!--- THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
	758	! PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
	759	! PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
	760	! 320 CONTINUE
	761
	762	!-----COMPUTE THE FRICTIONAL VELOCITY:
	763	! ZA(1982) EQS(2.60),(2.61).
	764
	765	DO 330 I=its,ite
	766	DTG=THX(I,KL)-THGB(I)
	767	PSIX=GZ1OZ0(I)-PSIM(I)
	768	IF((XLAND(I)-1.5).GE.0)THEN
	769	ZL=ZNT(I)
	770	ELSE
	771	ZL=0.01
	772	ENDIF
	773	PSIQ=ALOG(KARMANUST(I)ZA(I,KL)/XKA+ZA(I,KL)/ZL)-PSIH(I)
	774	UST(I)=KARMAN*WSPD(I)/PSIX
	775	!
	776	USTM=AMAX1(UST(I),0.1)
	777	IF((XLAND(I)-1.5).GE.0)THEN
	778	UST(I)=UST(I)
	779	ELSE
	780	UST(I)=USTM
	781	ENDIF
	782	! MOL(I,J)=KARMAN*DTG/(GZ1OZ0(I)-PSIH(I))/PRT
	783	330 CONTINUE
	784	!
	785	DO 420 I=its,ite
	786	WSPD1(I)=SQRT(UX(I,KL)UX(I,KL)+VX(I,KL)VX(I,KL))+1.E-9
	787	420 CONTINUE
	788	!
	789	!---- COMPUTE VERTICAL DIFFUSION
	790	!
	791	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	792	! COMPUTE PRELIMINARY VARIABLES
	793	!
	794	!
	795	DT4=2.*DTSTEP
	796	RDT=1./DT4
	797	!
	798	DO I=its,ite
	799	HGAMT(I)=0.
	800	HGAMQ(I)=0.
	801	WSCALE(I)=0.
	802	KPBL(I)=KL
	803	PBL(I)=ZQ(I,KL)
	804	KPBL0(I)=KL
	805	PBL0(I)=ZQ(I,KL)
	806	PBLFLG(I)=.TRUE.
	807	SFCFLG(I)=.TRUE.
	808	IF(BR(I).GT.0.0)SFCFLG(I)=.FALSE.
	809	ZL1(I)=ZA(I,KL)
	810	THERMAL(I)=THVX(I,KL)
	811	ENDDO
	812
	813	! COMPUTE THE FIRST GUESS OF PBL HEIGHT
	814
	815	DO I=its,ite
	816	STABLE(I)=.FALSE.
	817	BRUP(I)=BR(I)
	818	ENDDO
	819	DO K=KLM,KLPBL,-1
	820	DO I=its,ite
	821	IF(.NOT.STABLE(I))THEN
	822	BRDN(I)=BRUP(I)
	823	SPDK2=MAX(UX(I,K)2+VX(I,K)2,1.)
	824	BRUP(I)=(THVX(I,K)-THERMAL(I))(GZA(I,K)/THVX(I,KL))/SPDK2
	825	KPBL(I)=K
	826	STABLE(I)=BRUP(I).GT.BRCR
	827	ENDIF
	828	ENDDO
	829	ENDDO
	830	!
	831	DO I=its,ite
	832	K=KPBL(I)
	833	IF(BRDN(I).GE.BRCR)THEN
	834	BRINT=0.
	835	ELSEIF(BRUP(I).LE.BRCR)THEN
	836	BRINT=1.
	837	ELSE
	838	BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))
	839	ENDIF
	840	PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
	841	IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1
	842	ENDDO
	843	!
	844	DO I=its,ite
	845	FM=GZ1OZ0(I)-PSIM(I)
	846	FH=GZ1OZ0(I)-PSIH(I)
	847	HOL(I)=MAX(BR(I)FMFM/FH,RIMIN)
	848	IF(SFCFLG(I))THEN
	849	HOL(I)=MIN(HOL(I),-ZFMIN)
	850	ELSE
	851	HOL(I)=MAX(HOL(I),ZFMIN)
	852	ENDIF
	853	!
	854	HOL1=HOL(I)PBL(I)/ZL1(I)SFCFRAC
	855	HOL(I)=-HOL(I)*PBL(I)/ZL1(I)
	856	IF(SFCFLG(I))THEN
	857	PHIM(I)=(1.-APHI16HOL1)*(-1./4.)
	858	PHIH(I)=(1.-APHI16HOL1)*(-1./2.)
	859	ELSE
	860	PHIM(I)=(1.+APHI5*HOL1)
	861	PHIH(I)=PHIM(I)
	862	ENDIF
	863	WSCALE(I)=UST(I)/PHIM(I)
	864	WSCALE(I)=MIN(WSCALE(I),UST(I)*APHI16)
	865	WSCALE(I)=MAX(WSCALE(I),UST(I)/APHI5)
	866	ENDDO
	867
	868	! COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION
	869	! UNDER UNSTABLE CONDITIONS
	870
	871	DO I=its,ite
	872	IF(SFCFLG(I))THEN
	873	GAMFAC=CFAC/RHOX(I)/WSCALE(I)
	874	HGAMT(I)=MIN(GAMFAC*HFX(I)/CPM(I),GAMCRT)
	875	HGAMQ(I)=MIN(GAMFAC*QFX(I),GAMCRQ)
	876	IF((XLAND(I)-1.5).GE.0)HGAMQ(I)=0.
	877	VPERT=HGAMT(I)+EP1THX(I,KL)HGAMQ(I)
	878	VPERT=MIN(VPERT,GAMCRT)
	879	THERMAL(I)=THERMAL(I)+MAX(VPERT,0.)
	880	HGAMT(I)=MAX(HGAMT(I),0.0)
	881	HGAMQ(I)=MAX(HGAMQ(I),0.0)
	882	ELSE
	883	PBLFLG(I)=.FALSE.
	884	ENDIF
	885	ENDDO
	886	!
	887	DO I=its,ite
	888	IF(PBLFLG(I))THEN
	889	KPBL(I)=KL
	890	PBL(I)=ZQ(I,KL)
	891	ENDIF
	892	ENDDO
	893	!
	894	! ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL
	895	!
	896	DO I=its,ite
	897	IF(PBLFLG(I))THEN
	898	STABLE(I)=.FALSE.
	899	BRUP(I)=BR(I)
	900	ENDIF
	901	ENDDO
	902	DO K=KLM,KLPBL,-1
	903	DO I=its,ite
	904	IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN
	905	BRDN(I)=BRUP(I)
	906	SPDK2=MAX((UX(I,K)2+VX(I,K)2),1.)
	907	BRUP(I)=(THVX(I,K)-THERMAL(I))(GZA(I,K)/THVX(I,KL))/SPDK2
	908	KPBL(I)=K
	909	STABLE(I)=BRUP(I).GT.BRCR
	910	ENDIF
	911	ENDDO
	912	ENDDO
	913	!
	914	DO I=its,ite
	915	IF(PBLFLG(I))THEN
	916	K=KPBL(I)
	917	IF(BRDN(I).GE.BRCR)THEN
	918	BRINT=0.
	919	ELSEIF(BRUP(I).LE.BRCR)THEN
	920	BRINT=1.
	921	ELSE
	922	BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))
	923	ENDIF
	924	PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
	925	IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1
	926	IF(KPBL(I).LE.1)PBLFLG(I)=.FALSE.
	927	ENDIF
	928	ENDDO
	929	!
	930	! DIAGNOSTIC PBL HEIGHT WITH BRCR EFFECTIVELY ZERO (PBL0)
	931	!
	932	DO I=its,ite
	933	IF(PBLFLG(I))THEN
	934	STABLE(I)=.FALSE.
	935	BRUP(I)=BR(I)
	936	ENDIF
	937	ENDDO
	938	DO K=KLM,KLPBL,-1
	939	DO I=its,ite
	940	IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN
	941	BRDN(I)=BRUP(I)
	942	SPDK2=MAX((UX(I,K)2+VX(I,K)2),1.)
	943	BRUP(I)=(THVX(I,K)-THERMAL(I))(GZA(I,K)/THVX(I,KL))/SPDK2
	944	KPBL0(I)=K
	945	STABLE(I)=BRUP(I).GT.0.0
	946	ENDIF
	947
	948	ENDDO
	949	ENDDO
	950	!
	951	DO I=its,ite
	952	IF(PBLFLG(I))THEN
	953	K=KPBL0(I)
	954	IF(BRDN(I).GE.0.0)THEN
	955	BRINT=0.
	956	ELSEIF(BRUP(I).LE.0.0)THEN
	957	BRINT=1.
	958	ELSE
	959	BRINT=(0.0-BRDN(I))/(BRUP(I)-BRDN(I))
	960	ENDIF
	961	PBL0(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
	962	IF(PBL0(I).LT.ZQ(I,KPBL0(I)+1))KPBL0(I)=KPBL0(I)+1
	963	IF(KPBL0(I).LE.1)PBLFLG(I)=.FALSE.
	964	ENDIF
	965	ENDDO
	966
	967	!
	968	! COMPUTE DIFFUSION COEFFICIENTS BELOW PBL
	969	!
	970	DO K=kte,KLPBL,-1
	971	DO I=its,ite
	972	IF(KPBL(I).LT.K)THEN
	973	PRNUM=(PHIH(I)/PHIM(I)+CFACKARMANSFCFRAC)
	974	PRNUM=MIN(PRNUM,PRMAX)
	975	PRNUM=MAX(PRNUM,PRMIN)
	976	ZFAC=MAX((1.-(ZQ(I,K)-ZL1(I))/(PBL(I)-ZL1(I))),ZFMIN)
	977	XKZO=CKZ*DZA(I,K-1)
	978	XKZM(I,K)=XKZO+WSCALE(I)KARMANZQ(I,K)ZFAC*PFAC
	979	XKZH(I,K)=XKZM(I,K)/PRNUM
	980	XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)
	981	XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)
	982	XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)
	983	XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)
	984	ENDIF
	985	ENDDO
	986	ENDDO
	987	!
	988	! COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE)
	989	!
	990	DO K=kts+1,kte
	991	DO I=its,ite
	992	XKZO=CKZ*DZA(I,K-1)
	993	IF(K.LE.KPBL(I))THEN
	994	SS=((UX(I,K-1)-UX(I,K))*(UX(I,K-1)-UX(I,K))+(VX(I,K-1)- &
	995	VX(I,K))(VX(I,K-1)-VX(I,K)))/(DZA(I,K-1)DZA(I,K-1))+ &
	996	1.E-9
	997	RI=GOVRTH(I)(THVX(I,K-1)-THVX(I,K))/(SSDZA(I,K-1))
	998	IF(IMVDIF.EQ.1)THEN
	999	IF((QCX(I,K)+QIX(I,K)).GT.0.01E-3.AND.(QCX(I,K-1)+ &
	1000	QIX(I,K-1)).GT.0.01E-3)THEN
	1001	! IN CLOUD
	1002	QMEAN=0.5*(QX(I,K)+QX(I,K-1))
	1003	TMEAN=0.5*(SCR3(I,K)+SCR3(I,K-1))
	1004	ALPH=XLV*QMEAN/R/TMEAN
	1005	CHI=XLVXLVQMEAN/CP/RV/TMEAN/TMEAN
	1006	RI=(1.+ALPH)(RI-GG/SS/TMEAN/CP*((CHI-ALPH)/(1.+CHI)))
	1007	ENDIF
	1008	ENDIF
	1009	ZK=KARMAN*ZQ(I,K)
	1010	RL2=(ZKRLAM/(RLAM+ZK))*2
	1011	DK=RL2*SQRT(SS)
	1012	IF(RI.LT.0.)THEN
	1013	! UNSTABLE REGIME
	1014	SRI=SQRT(-RI)
	1015	XKZM(I,K)=XKZO+DK(1+8.(-RI)/(1+1.746*SRI))
	1016	XKZH(I,K)=XKZO+DK(1+8.(-RI)/(1+1.286*SRI))
	1017	ELSE
	1018	! STABLE REGIME
	1019	XKZH(I,K)=XKZO+DK/(1+5.RI)*2
	1020	PRNUM=1.0+2.1*RI
	1021	PRNUM=MIN(PRNUM,PRMAX)
	1022	XKZM(I,K)=(XKZH(I,K)-XKZO)*PRNUM+XKZO
	1023	ENDIF
	1024	!
	1025	XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)
	1026	XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)
	1027	XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)
	1028	XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)
	1029	ENDIF
	1030	!
	1031	ENDDO
	1032	ENDDO
	1033
	1034	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE
	1035
	1036	DO I=its,ite
	1037	DO K=kts,kte
	1038	AU(i,k)=0.
	1039	AL(i,k)=0.
	1040	AD(i,k)=0.
	1041	A1(i,k)=0.
	1042	A2(i,k)=0.
	1043	ENDDO
	1044	ENDDO
	1045
	1046	DO I=its,ite
	1047	AD(I,1)=1.
	1048	A1(I,1)=SCR3(I,KL)+HFX(I)/(RHOX(I)CPM(I))/ZQ(I,KL)DT4
	1049	A2(I,1)=QX(I,KL)+QFX(I)/(RHOX(I))/ZQ(I,KL)*DT4
	1050	ENDDO
	1051	!
	1052	DO K=kte,kts+1,-1
	1053	KK=kme-K
	1054	DO I=its,ite
	1055	DTODSD=DT4/dz8w2d(I,K)
	1056	DTODSU=DT4/dz8w2d(I,K-1)
	1057	DSIG=z2d(I,K)-z2d(I,K-1)
	1058	DSIG=-DSIG
	1059	RDZ=1./DZA(I,K-1)
	1060	IF(PBLFLG(I).AND.KPBL(I).LT.K)THEN
	1061	DSDZT=DSIGXKZH(I,K)RDZ*(G/CP-HGAMT(I)/PBL(I))
	1062	DSDZQ=DSIGXKZH(I,K)RDZ*(-HGAMQ(I)/PBL(I))
	1063	A2(I,KK)=A2(I,KK)+DTODSD*DSDZQ
	1064	A2(I,KK+1)=QX(I,K-1)-DTODSU*DSDZQ
	1065	ELSE
	1066	DSDZT=DSIGXKZH(I,K)RDZ*(G/CP)
	1067	A2(I,KK+1)=QX(I,K-1)
	1068	ENDIF
	1069	DSDZ2=DSIGXKZH(I,K)RDZ*RDZ
	1070	AU(I,KK)=-DTODSD*DSDZ2
	1071	AL(I,KK)=-DTODSU*DSDZ2
	1072	AD(I,KK)=AD(I,KK)-AU(I,KK)
	1073	AD(I,KK+1)=1.-AL(I,KK)
	1074	A1(I,KK)=A1(I,KK)+DTODSD*DSDZT
	1075	A1(I,KK+1)=SCR3(I,K-1)-DTODSU*DSDZT
	1076	ENDDO
	1077	ENDDO
	1078
	1079	! SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE
	1080
	1081	CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
	1082	its,ite,kts,kte )
	1083
	1084	! RECOVER TENDENCIES OF HEAT AND MOISTURE
	1085
	1086	DO K=kte,kts,-1
	1087	KK=kme-K
	1088	DO I=its,ite
	1089	TTEND=(A1(I,KK)-SCR3(I,K))*RDT
	1090	QTEND=(A2(I,KK)-QX(I,K))*RDT
	1091	TTNP(I,K)=TTNP(I,K)+TTEND
	1092	QTNP(I,K)=QTNP(I,K)+QTEND
	1093	ENDDO
	1094	ENDDO
	1095
	1096	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM
	1097
	1098	DO I=its,ite
	1099	DO K=kts,kte
	1100	AU(i,k)=0.
	1101	AL(i,k)=0.
	1102	AD(i,k)=0.
	1103	A1(i,k)=0.
	1104	A2(i,k)=0.
	1105	ENDDO
	1106	ENDDO
	1107
	1108	DO I=its,ite
	1109	AD(I,1)=1.
	1110	A1(I,1)=UX(I,KL)-UX(I,KL)/WSPD1(I)UST(I)UST(I)/ZQ(I,KL) &
	1111	DT4(WSPD1(I)/WSPD(I))**2
	1112	A2(I,1)=VX(I,KL)-VX(I,KL)/WSPD1(I)UST(I)UST(I)/ZQ(I,KL) &
	1113	DT4(WSPD1(I)/WSPD(I))**2
	1114	ENDDO
	1115	!
	1116	DO K=kte,kts+1,-1
	1117	KK=kme-K
	1118	DO I=its,ite
	1119	DTODSD=DT4/dz8w2d(I,K)
	1120	DTODSU=DT4/dz8w2d(I,K-1)
	1121	DSIG=z2d(I,K)-z2d(I,K-1)
	1122	DSIG=-DSIG
	1123	RDZ=1./DZA(I,K-1)
	1124	DSDZ2=DSIGXKZM(I,K)RDZ*RDZ
	1125	AU(I,KK)=-DTODSD*DSDZ2
	1126	AL(I,KK)=-DTODSU*DSDZ2
	1127	AD(I,KK)=AD(I,KK)-AU(I,KK)
	1128	AD(I,KK+1)=1.-AL(I,KK)
	1129	A1(I,KK+1)=UX(I,K-1)
	1130	A2(I,KK+1)=VX(I,K-1)
	1131	ENDDO
	1132	ENDDO
	1133
	1134	! SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM
	1135
	1136	CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
	1137	its,ite,kts,kte )
	1138
	1139	! RECOVER TENDENCIES OF MOMENTUM
	1140
	1141	DO K=kte,kts,-1
	1142	KK=kme-K
	1143	DO I=its,ite
	1144	UTEND=(A1(I,KK)-UX(I,K))*RDT
	1145	VTEND=(A2(I,KK)-VX(I,K))*RDT
	1146	UTNP(I,K)=UTNP(I,K)+UTEND
	1147	VTNP(I,K)=VTNP(I,K)+VTEND
	1148	ENDDO
	1149	ENDDO
	1150
	1151	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR CLOUD
	1152
	1153	DO I=its,ite
	1154	DO K=kts,kte
	1155	AU(i,k)=0.
	1156	AL(i,k)=0.
	1157	AD(i,k)=0.
	1158	A1(i,k)=0.
	1159	A2(i,k)=0.
	1160	ENDDO
	1161	ENDDO
	1162
	1163	! IF(IMOIST.EQ.1)GOTO 690
	1164	DO I=its,ite
	1165	AD(I,1)=1.
	1166	A1(I,1)=QCX(I,KL)
	1167	A2(I,1)=QIX(I,KL)
	1168	ENDDO
	1169	!
	1170	DO K=kte,kts+1,-1
	1171	KK=kme-K
	1172	DO I=its,ite
	1173	DTODSD=DT4/dz8w2d(I,K)
	1174	DTODSU=DT4/dz8w2d(I,K-1)
	1175	DSIG=z2d(I,K)-z2d(I,K-1)
	1176	DSIG=-DSIG
	1177	RDZ=1./DZA(I,K-1)
	1178	A1(I,KK+1)=QCX(I,K-1)
	1179	A2(I,KK+1)=QIX(I,K-1)
	1180	DSDZ2=DSIGXKZH(I,K)RDZ*RDZ
	1181	AU(I,KK)=-DTODSD*DSDZ2
	1182	AL(I,KK)=-DTODSU*DSDZ2
	1183	AD(I,KK)=AD(I,KK)-AU(I,KK)
	1184	AD(I,KK+1)=1.-AL(I,KK)
	1185	ENDDO
	1186	ENDDO
	1187
	1188	! SOLVE TRIDIAGONAL PROBLEM FOR CLOUD
	1189
	1190	CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
	1191	its,ite,kts,kte )
	1192	!
	1193	DO K=kte,kts,-1
	1194	KK=kme-K
	1195	DO I=its,ite
	1196	QCTEND=(A1(I,KK)-QCX(I,K))*RDT
	1197	QITEND=(A2(I,KK)-QIX(I,K))*RDT
	1198	QCTNP(I,K)=QCTNP(I,K)+QCTEND
	1199	QITNP(I,K)=QITNP(I,K)+QITEND
	1200	ENDDO
	1201	ENDDO
	1202	!
	1203	!---- END OF VERTICAL DIFFUSION
	1204	!
	1205	690 CONTINUE
	1206	!
	1207	!-----CALCULATION OF NEW VALUES DUE TO VERTICAL EXCHANGE PROCESSES IS
	1208	! COMPLETED. THE FINAL STEP IS TO ADD THE TENDENCIES CALCULATED
	1209	! IN HIRPBL TO THOSE OF MM4.
	1210
	1211	DO 820 K=kts,kte
	1212	NK=kme-K
	1213	DO 820 I=its,ite
	1214	U2DTEN(I,NK)=UTNP(I,K)
	1215	V2DTEN(I,NK)=VTNP(I,K)
	1216	820 CONTINUE
	1217	!
	1218	! IF(J.EQ.1.AND.IN.GT.1)GOTO 860
	1219	!SUE JBGN=3
	1220	!SUE JEND=JLXM-1
	1221
	1222	! change when nest
	1223	!SUE JBGN=2
	1224	!SUE JEND=JLXM
	1225
	1226	JBGN=jts
	1227	JEND=jte
	1228	IBGN=its
	1229	IEND=ite
	1230
	1231	! IF(J.LT.JBGN.OR.J.GT.JEND)GOTO 860
	1232	!SUE IBGN=3
	1233	!SUE IEND=ILXM-1
	1234
	1235	! change when nest
	1236	!SUE IBGN=2
	1237	!SUE IEND=ILXM
	1238
	1239	DO 830 K=kts,kte
	1240	NK=kme-K
	1241	DO 830 I=IBGN,IEND
	1242	T2DTEN(I,NK)=TTNP(I,K)
	1243	830 CONTINUE
	1244	!
	1245	! IF(IDRY.EQ.1)GOTO 860
	1246	DO 840 K=kts,kte
	1247	NK=kme-K
	1248	DO 840 I=IBGN,IEND
	1249	QV2DTEN(I,NK)=QTNP(I,K)
	1250	840 CONTINUE
	1251
	1252	! IF(IMOIST.EQ.1)GOTO 860
	1253	DO 850 K=kts,kte
	1254	NK=kme-K
	1255	DO 850 I=IBGN,IEND
	1256	QC2DTEN(I,NK)=QCTNP(I,K)
	1257	850 CONTINUE
	1258
	1259	IF(flag_QI .AND. PRESENT( QI2DTEN ) ) THEN
	1260	DO K=kts,kte
	1261	NK=kme-K
	1262	DO I=IBGN,IEND
	1263	QI2DTEN(I,NK)=QITNP(I,K)
	1264	ENDDO
	1265	ENDDO
	1266	ENDIF
	1267
	1268	860 CONTINUE
	1269	!
	1270	!-----APPLY ASSELIN FILTER TO TGD FOR LARGE TIME STEP:
	1271	!
	1272	! DO 885 I=its,ite
	1273	! TSK(I)=TSK(I)(PS(I)/100.)*ROVCP
	1274	! 885 CONTINUE
	1275	!
	1276	940 CONTINUE
	1277	!
	1278	! KPBL IS NEEDED FOR THE FDDA, AND SINCE THERE IS NO LONGER JUST ONE
	1279	! LARGE "J LOOP" IT MUST BE STORED AS (I,J)...
	1280	!
	1281	! USE NEW DIAGNOSED PBL DEPTH (CALCULATED WITH brcr=0.0)
	1282	! PBL IS USED FOR OUTPUT AND NEXT-TIME-STEP BELJAARS CALC IN SFCLAY
	1283	DO 950 I=its,ite
	1284	KPBL1D(I)=KPBL0(I)
	1285	PBL(I)=PBL0(I)
	1286	950 CONTINUE
	1287
	1288	END SUBROUTINE MRF2D
	1289
	1290	!================================================================
	1291	SUBROUTINE TRIDI2(CL,CM,CU,R1,R2,AU,A1,A2, &
	1292	its,ite,kts,kte )
	1293	!----------------------------------------------------------------
	1294	IMPLICIT NONE
	1295	!----------------------------------------------------------------
	1296
	1297	INTEGER, INTENT(IN ) :: its,ite, kts,kte
	1298
	1299	REAL, DIMENSION( its:ite, kts+1:kte+1 ) , &
	1300	INTENT(IN ) :: CL
	1301
	1302	REAL, DIMENSION( its:ite, kts:kte ) , &
	1303	INTENT(IN ) :: CM, &
	1304	R1, &
	1305	R2
	1306	REAL, DIMENSION( its:ite, kts:kte ) , &
	1307	INTENT(INOUT) :: AU, &
	1308	CU, &
	1309	A1, &
	1310	A2
	1311
	1312	REAL :: FK
	1313	INTEGER :: I,K,L,N
	1314
	1315	!----------------------------------------------------------------
	1316
	1317	L=ite
	1318	N=kte
	1319
	1320	DO I=its,L
	1321	FK=1./CM(I,1)
	1322	AU(I,1)=FK*CU(I,1)
	1323	A1(I,1)=FK*R1(I,1)
	1324	A2(I,1)=FK*R2(I,1)
	1325	ENDDO
	1326	DO K=2,N-1
	1327	DO I=its,L
	1328	FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
	1329	AU(I,K)=FK*CU(I,K)
	1330	A1(I,K)=FK(R1(I,K)-CL(I,K)A1(I,K-1))
	1331	A2(I,K)=FK(R2(I,K)-CL(I,K)A2(I,K-1))
	1332	ENDDO
	1333	ENDDO
	1334	DO I=its,L
	1335	FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
	1336	A1(I,N)=FK(R1(I,N)-CL(I,N)A1(I,N-1))
	1337	A2(I,N)=FK(R2(I,N)-CL(I,N)A2(I,N-1))
	1338
	1339	ENDDO
	1340	DO K=N-1,kts,-1
	1341	DO I=its,L
	1342	A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
	1343	A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
	1344	ENDDO
	1345	ENDDO
	1346
	1347	END SUBROUTINE TRIDI2
	1348
	1349	!===================================================================
	1350	SUBROUTINE mrfinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
	1351	RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &
	1352	restart, allowed_to_read , &
	1353	ids, ide, jds, jde, kds, kde, &
	1354	ims, ime, jms, jme, kms, kme, &
	1355	its, ite, jts, jte, kts, kte )
	1356	!-------------------------------------------------------------------
	1357	IMPLICIT NONE
	1358	!-------------------------------------------------------------------
	1359	LOGICAL , INTENT(IN) :: restart , allowed_to_read
	1360	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
	1361	ims, ime, jms, jme, kms, kme, &
	1362	its, ite, jts, jte, kts, kte
	1363	INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR
	1364
	1365	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
	1366	RUBLTEN, &
	1367	RVBLTEN, &
	1368	RTHBLTEN, &
	1369	RQVBLTEN, &
	1370	RQCBLTEN, &
	1371	RQIBLTEN
	1372	INTEGER :: i, j, k, itf, jtf, ktf
	1373
	1374	jtf=min0(jte,jde-1)
	1375	ktf=min0(kte,kde-1)
	1376	itf=min0(ite,ide-1)
	1377
	1378	IF(.not.restart)THEN
	1379	DO j=jts,jtf
	1380	DO k=kts,ktf
	1381	DO i=its,itf
	1382	RUBLTEN(i,k,j)=0.
	1383	RVBLTEN(i,k,j)=0.
	1384	RTHBLTEN(i,k,j)=0.
	1385	RQVBLTEN(i,k,j)=0.
	1386	RQCBLTEN(i,k,j)=0.
	1387	ENDDO
	1388	ENDDO
	1389	ENDDO
	1390	ENDIF
	1391
	1392	IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
	1393	DO j=jts,jtf
	1394	DO k=kts,ktf
	1395	DO i=its,itf
	1396	RQIBLTEN(i,k,j)=0.
	1397	ENDDO
	1398	ENDDO
	1399	ENDDO
	1400	ENDIF
	1401
	1402	END SUBROUTINE mrfinit
	1403
	1404	!-------------------------------------------------------------------
	1405
	1406	END MODULE module_bl_mrf
	1407

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

Download in other formats:

Original Format