Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_bl_gfs.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: 40.7 KB

Rev	Line
[2759]	1	!LWRF:MODEL_LAYER:PHYSICS
	2	!
	3	MODULE module_bl_gfs
	4
	5	CONTAINS
	6
	7	!-------------------------------------------------------------------
	8	SUBROUTINE BL_GFS(U3D,V3D,TH3D,T3D,QV3D,QC3D, P3D,PI3D, &
	9	RUBLTEN,RVBLTEN,RTHBLTEN, &
	10	RQVBLTEN,RQCBLTEN, &
	11	CP,G,ROVCP,R,ROVG,FLAG_QI, &
	12	dz8w,z,PSFC, &
	13	UST,PBL,PSIM,PSIH, &
	14	HFX,QFX,TSK,GZ1OZ0,WSPD,BR, &
	15	DT,KPBL2D,EP1,KARMAN, &
	16	ids,ide, jds,jde, kds,kde, &
	17	ims,ime, jms,jme, kms,kme, &
	18	its,ite, jts,jte, kts,kte, &
	19	! optional
	20	qi3d,rqiblten )
	21	!--------------------------------------------------------------------
	22	USE MODULE_GFS_MACHINE, ONLY : kind_phys
	23	!-------------------------------------------------------------------
	24	IMPLICIT NONE
	25	!-------------------------------------------------------------------
	26	!-- U3D 3D u-velocity interpolated to theta points (m/s)
	27	!-- V3D 3D v-velocity interpolated to theta points (m/s)
	28	!-- TH3D 3D potential temperature (K)
	29	!-- T3D temperature (K)
	30	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
	31	!-- QC3D 3D cloud mixing ratio (Kg/Kg)
	32	!-- QI3D 3D ice mixing ratio (Kg/Kg)
	33	!-- P3D 3D pressure (Pa)
	34	!-- PI3D 3D exner function (dimensionless)
	35	!-- rr3D 3D dry air density (kg/m^3)
	36	!-- RUBLTEN U tendency due to
	37	! PBL parameterization (m/s^2)
	38	!-- RVBLTEN V tendency due to
	39	! PBL parameterization (m/s^2)
	40	!-- RTHBLTEN Theta tendency due to
	41	! PBL parameterization (K/s)
	42	!-- RQVBLTEN Qv tendency due to
	43	! PBL parameterization (kg/kg/s)
	44	!-- RQCBLTEN Qc tendency due to
	45	! PBL parameterization (kg/kg/s)
	46	!-- RQIBLTEN Qi tendency due to
	47	! PBL parameterization (kg/kg/s)
	48	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
	49	!-- G acceleration due to gravity (m/s^2)
	50	!-- ROVCP R/CP
	51	!-- R gas constant for dry air (J/kg/K)
	52	!-- ROVG R/G
	53	!-- P_QI species index for cloud ice
	54	!-- dz8w dz between full levels (m)
	55	!-- z height above sea level (m)
	56	!-- PSFC pressure at the surface (Pa)
	57	!-- UST u* in similarity theory (m/s)
	58	!-- PBL PBL height (m)
	59	!-- PSIM similarity stability function for momentum
	60	!-- PSIH similarity stability function for heat
	61	!-- HFX upward heat flux at the surface (W/m^2)
	62	!-- QFX upward moisture flux at the surface (kg/m^2/s)
	63	!-- TSK surface temperature (K)
	64	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
	65	!-- WSPD wind speed at lowest model level (m/s)
	66	!-- BR bulk Richardson number in surface layer
	67	!-- DT time step (s)
	68	!-- rvovrd R_v divided by R_d (dimensionless)
	69	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
	70	!-- KARMAN Von Karman constant
	71	!-- ids start index for i in domain
	72	!-- ide end index for i in domain
	73	!-- jds start index for j in domain
	74	!-- jde end index for j in domain
	75	!-- kds start index for k in domain
	76	!-- kde end index for k in domain
	77	!-- ims start index for i in memory
	78	!-- ime end index for i in memory
	79	!-- jms start index for j in memory
	80	!-- jme end index for j in memory
	81	!-- kms start index for k in memory
	82	!-- kme end index for k in memory
	83	!-- its start index for i in tile
	84	!-- ite end index for i in tile
	85	!-- jts start index for j in tile
	86	!-- jte end index for j in tile
	87	!-- kts start index for k in tile
	88	!-- kte end index for k in tile
	89	!-------------------------------------------------------------------
	90
	91	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
	92	ims,ime, jms,jme, kms,kme, &
	93	its,ite, jts,jte, kts,kte
	94
	95	LOGICAL, INTENT(IN) :: flag_qi
	96
	97	REAL, INTENT(IN) :: &
	98	CP, &
	99	DT, &
	100	EP1, &
	101	G, &
	102	KARMAN, &
	103	R, &
	104	ROVCP, &
	105	ROVG
	106
	107	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
	108	DZ8W, &
	109	P3D, &
	110	PI3D, &
	111	QC3D, &
	112	QV3D, &
	113	T3D, &
	114	TH3D, &
	115	U3D, &
	116	V3D, &
	117	Z
	118
	119
	120	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: &
	121	RTHBLTEN, &
	122	RQCBLTEN, &
	123	RQVBLTEN, &
	124	RUBLTEN, &
	125	RVBLTEN
	126
	127	REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
	128	BR, &
	129	GZ1OZ0, &
	130	HFX, &
	131	PSFC, &
	132	PSIM, &
	133	PSIH, &
	134	QFX, &
	135	TSK
	136
	137	REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
	138	PBL, &
	139	UST, &
	140	WSPD
	141
	142	INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
	143	KPBL2D
	144
	145	!--------------------------- OPTIONAL VARS ------------------------------
	146	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), &
	147	OPTIONAL, INTENT(IN) :: &
	148	QI3D
	149
	150	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), &
	151	OPTIONAL, INTENT(INOUT) :: &
	152	RQIBLTEN
	153
	154	!--------------------------- LOCAL VARS ------------------------------
	155
	156
	157	REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
	158	DEL, &
	159	DU, &
	160	DV, &
	161	PHIL, &
	162	PRSL, &
	163	PRSLK, &
	164	T1, &
	165	TAU, &
	166	U1, &
	167	V1
	168
	169	REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
	170	PHII, &
	171	PRSI
	172
	173	REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 3) :: &
	174	Q1, &
	175	RTG
	176
	177	REAL (kind=kind_phys), DIMENSION(its:ite) :: &
	178	DQSFC, &
	179	DTSFC, &
	180	DUSFC, &
	181	DVSFC, &
	182	EVAP, &
	183	FH, &
	184	FM, &
	185	HEAT, &
	186	HGAMQ, &
	187	HGAMT, &
	188	HPBL, &
	189	PSK, &
	190	QSS, &
	191	RBSOIL, &
	192	RCL, &
	193	SPD1, &
	194	STRESS, &
	195	TSEA
	196
	197	REAL (kind=kind_phys) :: &
	198	CPM, &
	199	DELTIM, &
	200	FMTMP, &
	201	RRHOX
	202
	203	INTEGER, DIMENSION( its:ite ) :: &
	204	KPBL
	205
	206	INTEGER :: &
	207	I, &
	208	IM, &
	209	J, &
	210	K, &
	211	KM, &
	212	KTEM, &
	213	KTEP, &
	214	KX, &
	215	L, &
	216	NTRAC
	217
	218	IM=ITE-ITS+1
	219	KX=KTE-KTS+1
	220	KTEM=KTE-1
	221	KTEP=KTE+1
	222	NTRAC=2
	223	DELTIM=DT
	224	IF (flag_qi) NTRAC=3
	225
	226
	227	DO J=jts,jte
	228
	229	DO i=its,ite
	230	RRHOX=(RT3D(I,KTS,J)(1.+EP1*QV3D(I,KTS,J)))/PSFC(I,J)
	231	CPM=CP(1.+0.8QV3D(i,kts,j))
	232	FMTMP=GZ1OZ0(i,j)-PSIM(i,j)
	233	PSK(i)=(PSFC(i,j).00001)*ROVCP
	234	FM(i)=FMTMP
	235	FH(i)=GZ1OZ0(i,j)-PSIH(i,j)
	236	TSEA(i)=TSK(i,j)
	237	QSS(i)=QV3D(i,kts,j) ! not used in moninp so set to qv3d for now
	238	HEAT(i)=HFX(i,j)/CPM*RRHOX
	239	EVAP(i)=QFX(i,j)*RRHOX
	240	STRESS(i)=KARMANKARMANWSPD(i,j)WSPD(i,j)/(FMTMPFMTMP)
	241	SPD1(i)=WSPD(i,j)
	242	PRSI(i,kts)=PSFC(i,j)*.001
	243	PHII(I,kts)=0.
	244	RCL(i)=1.
	245	RBSOIL(I)=BR(i,j)
	246	ENDDO
	247
	248	DO k=kts,kte
	249	DO i=its,ite
	250	DV(I,K) = 0.
	251	DU(I,K) = 0.
	252	TAU(I,K) = 0.
	253	U1(I,K) = U3D(i,k,j)
	254	V1(I,K) = V3D(i,k,j)
	255	T1(I,K) = T3D(i,k,j)
	256	Q1(I,K,1) = QV3D(i,k,j)/(1.+QV3D(i,k,j))
	257	Q1(I,K,2) = QC3D(i,k,j)/(1.+QC3D(i,k,j))
	258	PRSL(I,K)=P3D(i,k,j)*.001
	259	ENDDO
	260	ENDDO
	261
	262	DO k=kts,kte
	263	DO i=its,ite
	264	PRSLK(I,K)=(PRSL(i,k).01)*ROVCP
	265	ENDDO
	266	ENDDO
	267
	268	DO k=kts+1,kte
	269	km=k-1
	270	DO i=its,ite
	271	DEL(i,km)=PRSL(i,km)/ROVG*dz8w(i,km,j)/T3D(i,km,j)
	272	PRSI(i,k)=PRSI(i,km)-DEL(i,km)
	273	PHII(I,K)=(Z(i,k,j)-Z(i,kts,j))*G
	274	PHIL(I,KM)=0.5(Z(i,k,j)+Z(i,km,j)-2.Z(i,kts,j))*G
	275	ENDDO
	276	ENDDO
	277
	278	DO i=its,ite
	279	DEL(i,kte)=DEL(i,ktem)
	280	PRSI(i,ktep)=PRSI(i,kte)-DEL(i,ktem)
	281	PHII(I,KTEP)=PHII(I,KTE)+dz8w(i,kte,j)*G
	282	PHIL(I,KTE)=PHII(I,KTE)-PHIL(I,KTEM)+PHII(I,KTE)
	283	ENDDO
	284
	285	IF (flag_QI .AND. PRESENT( QI3D ) ) THEN
	286	DO k=kts,kte
	287	DO i=its,ite
	288	Q1(I,K,3) = QI3D(i,k,j)/(1.+QI3D(i,k,j))
	289	ENDDO
	290	ENDDO
	291	ENDIF
	292
	293	DO l=1,ntrac
	294	DO k=kts,kte
	295	DO i=its,ite
	296	RTG(I,K,L) = 0.
	297	ENDDO
	298	ENDDO
	299	ENDDO
	300
	301	CALL MONINP(IM,IM,KX,NTRAC,DV,DU,TAU,RTG,U1,V1,T1,Q1, &
	302	PSK,RBSOIL,FM,FH,TSEA,QSS,HEAT,EVAP,STRESS, &
	303	SPD1,KPBL,PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL, &
	304	DELTIM,DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,HGAMQ)
	305
	306
	307	DO k=kts,kte
	308	DO i=its,ite
	309	RVBLTEN(I,K,J)=DV(I,K)
	310	RUBLTEN(I,K,J)=DU(I,K)
	311	RTHBLTEN(I,K,J)=TAU(I,K)/PI3D(I,K,J)
	312	RQVBLTEN(I,K,J)=RTG(I,K,1)/(1.-Q1(I,K,1))**2
	313	RQCBLTEN(I,K,J)=RTG(I,K,2)/(1.-Q1(I,K,2))**2
	314	ENDDO
	315	ENDDO
	316
	317	IF (flag_QI .AND. PRESENT( RQIBLTEN )) THEN
	318	DO k=kts,kte
	319	DO i=its,ite
	320	RQIBLTEN(I,K,J)=RTG(I,K,3)/(1.-Q1(I,K,3))**2
	321	ENDDO
	322	ENDDO
	323	ENDIF
	324
	325	DO i=its,ite
	326	UST(i,j)=SQRT(STRESS(i))
	327	WSPD(i,j)=SQRT(U3D(I,KTS,J)*U3D(I,KTS,J)+ &
	328	V3D(I,KTS,J)*V3D(I,KTS,J))+1.E-9
	329	PBL(i,j)=HPBL(i)
	330	KPBL2D(i,j)=kpbl(i)
	331	ENDDO
	332
	333	ENDDO
	334
	335
	336	END SUBROUTINE BL_GFS
	337
	338	!===================================================================
	339	SUBROUTINE gfsinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
	340	RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &
	341	restart, &
	342	allowed_to_read, &
	343	ids, ide, jds, jde, kds, kde, &
	344	ims, ime, jms, jme, kms, kme, &
	345	its, ite, jts, jte, kts, kte )
	346	!-------------------------------------------------------------------
	347	IMPLICIT NONE
	348	!-------------------------------------------------------------------
	349	LOGICAL , INTENT(IN) :: allowed_to_read,restart
	350	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
	351	ims, ime, jms, jme, kms, kme, &
	352	its, ite, jts, jte, kts, kte
	353	INTEGER , INTENT(IN) :: P_QI, P_FIRST_SCALAR
	354
	355	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
	356	RUBLTEN, &
	357	RVBLTEN, &
	358	RTHBLTEN, &
	359	RQVBLTEN, &
	360	RQCBLTEN, &
	361	RQIBLTEN
	362	INTEGER :: i, j, k, itf, jtf, ktf
	363
	364	jtf=min0(jte,jde-1)
	365	ktf=min0(kte,kde-1)
	366	itf=min0(ite,ide-1)
	367
	368	IF(.not.restart)THEN
	369	DO j=jts,jtf
	370	DO k=kts,ktf
	371	DO i=its,itf
	372	RUBLTEN(i,k,j)=0.
	373	RVBLTEN(i,k,j)=0.
	374	RTHBLTEN(i,k,j)=0.
	375	RQVBLTEN(i,k,j)=0.
	376	RQCBLTEN(i,k,j)=0.
	377	ENDDO
	378	ENDDO
	379	ENDDO
	380	ENDIF
	381
	382	IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
	383	DO j=jts,jtf
	384	DO k=kts,ktf
	385	DO i=its,itf
	386	RQIBLTEN(i,k,j)=0.
	387	ENDDO
	388	ENDDO
	389	ENDDO
	390	ENDIF
	391
	392	END SUBROUTINE gfsinit
	393
	394	! --------------------------------------------------------------
	395	!FPP$ NOCONCUR R
	396	SUBROUTINE MONINP(IX,IM,KM,ntrac,DV,DU,TAU,RTG, &
	397	& U1,V1,T1,Q1, &
	398	& PSK,RBSOIL,FM,FH,TSEA,QSS,HEAT,EVAP,STRESS,SPD1,KPBL, &
	399	! & PSK,RBSOIL,CD,CH,FM,FH,TSEA,QSS,DPHI,SPD1,KPBL, &
	400	& PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,DELTIM, &
	401	& DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,HGAMQ)
	402	!
	403	USE MODULE_GFS_MACHINE, ONLY : kind_phys
	404	USE MODULE_GFS_PHYSCONS, grav => con_g, RD => con_RD, CP => con_CP &
	405	&, HVAP => con_HVAP, ROG => con_ROG, FV => con_FVirt
	406
	407	implicit none
	408	!
	409	! include 'constant.h'
	410	!
	411	!
	412	! Arguments
	413	!
	414	integer IX, IM, KM, ntrac, KPBL(IM)
	415	!
	416	real(kind=kind_phys) DELTIM
	417	real(kind=kind_phys) DV(IM,KM), DU(IM,KM), &
	418	& TAU(IM,KM), RTG(IM,KM,ntrac), &
	419	& U1(IX,KM), V1(IX,KM), &
	420	& T1(IX,KM), Q1(IX,KM,ntrac), &
	421	& PSK(IM), RBSOIL(IM), &
	422	! & CD(IM), CH(IM), &
	423	& FM(IM), FH(IM), &
	424	& TSEA(IM), QSS(IM), &
	425	& SPD1(IM), &
	426	! & DPHI(IM), SPD1(IM), &
	427	& PRSI(IX,KM+1), DEL(IX,KM), &
	428	& PRSL(IX,KM), PRSLK(IX,KM), &
	429	& PHII(IX,KM+1), PHIL(IX,KM), &
	430	& RCL(IM), DUSFC(IM), &
	431	& dvsfc(IM), dtsfc(IM), &
	432	& DQSFC(IM), HPBL(IM), &
	433	& HGAMT(IM), hgamq(IM)
	434	!
	435	! Locals
	436	!
	437	integer i,iprt,is,iun,k,kk,kmpbl,lond
	438	! real(kind=kind_phys) betaq(IM), betat(IM), betaw(IM), &
	439	real(kind=kind_phys) evap(IM), heat(IM), phih(IM), &
	440	& phim(IM), rbdn(IM), rbup(IM), &
	441	& the1(IM), stress(im), beta(im), &
	442	& the1v(IM), thekv(IM), thermal(IM), &
	443	& thesv(IM), ustar(IM), wscale(IM)
	444	! & thesv(IM), ustar(IM), wscale(IM), zl1(IM)
	445	!
	446	real(kind=kind_phys) RDZT(IM,KM-1), &
	447	& ZI(IM,KM+1), ZL(IM,KM), &
	448	& DKU(IM,KM-1), DKT(IM,KM-1), DKO(IM,KM-1), &
	449	& AL(IM,KM-1), AD(IM,KM), &
	450	& AU(IM,KM-1), A1(IM,KM), &
	451	& A2(IM,KM), THETA(IM,KM), &
	452	& AT(IM,KM*(ntrac-1))
	453	logical pblflg(IM), sfcflg(IM), stable(IM)
	454	!
	455	real(kind=kind_phys) aphi16, aphi5, bet1, bvf2, &
	456	& cfac, conq, cont, conw, &
	457	& conwrc, dk, dkmax, dkmin, &
	458	& dq1, dsdz2, dsdzq, dsdzt, &
	459	& dsig, dt, dthe1, dtodsd, &
	460	& dtodsu, dw2, dw2min, g, &
	461	& gamcrq, gamcrt, gocp, gor, gravi, &
	462	& hol, pfac, prmax, prmin, prinv, &
	463	& prnum, qmin, qtend, rbcr, &
	464	& rbint, rdt, rdz, rdzt1, &
	465	& ri, rimin, rl2, rlam, &
	466	& rone, rzero, sfcfrac, &
	467	& sflux, shr2, spdk2, sri, &
	468	& tem, ti, ttend, tvd, &
	469	& tvu, utend, vk, vk2, &
	470	& vpert, vtend, xkzo, zfac, &
	471	& zfmin, zk, tem1
	472	!
	473	PARAMETER(g=grav)
	474	PARAMETER(GOR=G/RD,GOCP=G/CP)
	475	PARAMETER(CONT=1000.CP/G,CONQ=1000.HVAP/G,CONW=1000./G)
	476	PARAMETER(RLAM=150.,VK=0.4,VK2=VK*VK,PRMIN=1.0,PRMAX=4.)
	477	PARAMETER(DW2MIN=0.0001,DKMIN=1.0,DKMAX=1000.,RIMIN=-100.)
	478	PARAMETER(RBCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1)
	479	PARAMETER(QMIN=1.E-8,XKZO=1.0,ZFMIN=1.E-8,APHI5=5.,APHI16=16.)
	480	! PARAMETER(GAMCRT=3.,GAMCRQ=2.E-3)
	481	PARAMETER(GAMCRT=3.,GAMCRQ=0.)
	482	PARAMETER(RZERO=0.,RONE=1.)
	483	PARAMETER(IUN=84)
	484	!
	485	!
	486	!-----------------------------------------------------------------------
	487	!
	488	601 FORMAT(1X,' MONINP LAT LON STEP HOUR ',3I6,F6.1)
	489	602 FORMAT(1X,' K',' Z',' T',' TH', &
	490	& ' TVH',' Q',' U',' V', &
	491	& ' SP')
	492	603 FORMAT(1X,I5,8F9.1)
	493	604 FORMAT(1X,' SFC',9X,F9.1,18X,F9.1)
	494	605 FORMAT(1X,' K ZL SPD2 THEKV THE1V' &
	495	& ,' THERMAL RBUP')
	496	606 FORMAT(1X,I5,6F8.2)
	497	607 FORMAT(1X,' KPBL HPBL FM FH HGAMT', &
	498	& ' HGAMQ WS USTAR CD CH')
	499	608 FORMAT(1X,I5,9F8.2)
	500	609 FORMAT(1X,' K PR DKT DKU ',I5,3F8.2)
	501	610 FORMAT(1X,' K PR DKT DKU ',I5,3F8.2,' L2 RI T2', &
	502	& ' SR2 ',2F8.2,2E10.2)
	503	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	504	! COMPUTE PRELIMINARY VARIABLES
	505	!
	506
	507	if (IX .lt. im) stop
	508	!
	509	IPRT = 0
	510	IF(IPRT.EQ.1) THEN
	511	!!! LATD = 0
	512	LOND = 0
	513	ELSE
	514	!!! LATD = 0
	515	LOND = 0
	516	ENDIF
	517	!
	518	gravi = 1.0 / grav
	519	DT = 2. * DELTIM
	520	RDT = 1. / DT
	521	KMPBL = KM / 2
	522	!
	523	do k=1,km
	524	do i=1,im
	525	zi(i,k) = phii(i,k) * gravi
	526	zl(i,k) = phil(i,k) * gravi
	527	enddo
	528	enddo
	529	!
	530	do k=1,kmpbl
	531	do i=1,im
	532	theta(i,k) = t1(i,k) * psk(i) / prslk(i,k)
	533	enddo
	534	enddo
	535	!
	536	DO K = 1,KM-1
	537	DO I=1,IM
	538	RDZT(I,K) = GOR * PRSI(I,K+1) / (PRSL(I,K) - PRSL(I,K+1))
	539	ENDDO
	540	ENDDO
	541	!
	542	DO I = 1,IM
	543	DUSFC(I) = 0.
	544	DVSFC(I) = 0.
	545	DTSFC(I) = 0.
	546	DQSFC(I) = 0.
	547	HGAMT(I) = 0.
	548	HGAMQ(I) = 0.
	549	WSCALE(I) = 0.
	550	KPBL(I) = 1
	551	HPBL(I) = ZI(I,2)
	552	PBLFLG(I) = .TRUE.
	553	SFCFLG(I) = .TRUE.
	554	IF(RBSOIL(I).GT.0.0) SFCFLG(I) = .FALSE.
	555	ENDDO
	556	!!
	557	DO I=1,IM
	558	RDZT1 = GOR * prSL(i,1) / DEL(i,1)
	559	! BET1 = DTRDZT1SPD1(I)/T1(I,1)
	560	BETA(I) = DT*RDZT1/T1(I,1)
	561	! BETAW(I) = BET1*CD(I)
	562	! BETAT(I) = BET1*CH(I)
	563	! BETAQ(I) = DPHI(I)*BETAT(I)
	564	ENDDO
	565	!
	566	DO I=1,IM
	567	! ZL1(i) = 0.-(T1(I,1)+TSEA(I))/2.LOG(PRSL(I,1)/PRSI(I,1))ROG
	568	! USTAR(I) = SQRT(CD(I)SPD1(I)*2)
	569	USTAR(I) = SQRT(STRESS(I))
	570	ENDDO
	571	!
	572	DO I=1,IM
	573	THESV(I) = TSEA(I)(1.+FVMAX(QSS(I),QMIN))
	574	THE1(I) = THETA(I,1)
	575	THE1V(I) = THE1(I)(1.+FVMAX(Q1(I,1,1),QMIN))
	576	THERMAL(I) = THE1V(I)
	577	! DTHE1 = (THE1(I)-TSEA(I))
	578	! DQ1 = (MAX(Q1(I,1,1),QMIN) - MAX(QSS(I),QMIN))
	579	! HEAT(I) = -CH(I)SPD1(I)DTHE1
	580	! EVAP(I) = -CH(I)SPD1(I)DQ1
	581	ENDDO
	582	!
	583	!
	584	! COMPUTE THE FIRST GUESS OF PBL HEIGHT
	585	!
	586	DO I=1,IM
	587	STABLE(I) = .FALSE.
	588	! ZL(i,1) = ZL1(i)
	589	RBUP(I) = RBSOIL(I)
	590	ENDDO
	591	DO K = 2, KMPBL
	592	DO I = 1, IM
	593	IF(.NOT.STABLE(I)) THEN
	594	RBDN(I) = RBUP(I)
	595	! ZL(I,k) = ZL(I,K-1) - (T1(i,k)+T1(i,K-1))/2 *
	596	! & LOG(PRSL(I,K)/PRSL(I,K-1)) * ROG
	597	THEKV(I) = THETA(i,k)(1.+FVMAX(Q1(i,k,1),QMIN))
	598	SPDK2 = MAX(RCL(i)(U1(i,k)U1(i,k)+V1(i,k)*V1(i,k)),RONE)
	599	RBUP(I) = (THEKV(I)-THE1V(I))(GZL(I,k)/THE1V(I))/SPDK2
	600	KPBL(I) = K
	601	STABLE(I) = RBUP(I).GT.RBCR
	602	ENDIF
	603	ENDDO
	604	ENDDO
	605	!
	606	DO I = 1,IM
	607	K = KPBL(I)
	608	IF(RBDN(I).GE.RBCR) THEN
	609	RBINT = 0.
	610	ELSEIF(RBUP(I).LE.RBCR) THEN
	611	RBINT = 1.
	612	ELSE
	613	RBINT = (RBCR-RBDN(I))/(RBUP(I)-RBDN(I))
	614	ENDIF
	615	HPBL(I) = ZL(I,K-1) + RBINT*(ZL(I,K)-ZL(I,K-1))
	616	IF(HPBL(I).LT.ZI(I,KPBL(I))) KPBL(I) = KPBL(I) - 1
	617	ENDDO
	618	!!
	619	DO I=1,IM
	620	HOL = MAX(RBSOIL(I)FM(I)FM(I)/FH(I),RIMIN)
	621	IF(SFCFLG(I)) THEN
	622	HOL = MIN(HOL,-ZFMIN)
	623	ELSE
	624	HOL = MAX(HOL,ZFMIN)
	625	ENDIF
	626	!
	627	! HOL = HOLHPBL(I)/ZL1(I)SFCFRAC
	628	HOL = HOLHPBL(I)/ZL(I,1)SFCFRAC
	629	IF(SFCFLG(I)) THEN
	630	! PHIM = (1.-APHI16HOL)*(-1./4.)
	631	! PHIH = (1.-APHI16HOL)*(-1./2.)
	632	TEM = 1.0 / (1. - APHI16*HOL)
	633	PHIH(I) = SQRT(TEM)
	634	PHIM(I) = SQRT(PHIH(I))
	635	ELSE
	636	PHIM(I) = (1.+APHI5*HOL)
	637	PHIH(I) = PHIM(I)
	638	ENDIF
	639	WSCALE(I) = USTAR(I)/PHIM(I)
	640	WSCALE(I) = MIN(WSCALE(I),USTAR(I)*APHI16)
	641	WSCALE(I) = MAX(WSCALE(I),USTAR(I)/APHI5)
	642	ENDDO
	643	!
	644	! COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION
	645	! UNDER UNSTABLE CONDITIONS
	646	!
	647	DO I = 1,IM
	648	SFLUX = HEAT(I) + EVAP(I)FVTHE1(I)
	649	IF(SFCFLG(I).AND.SFLUX.GT.0.0) THEN
	650	HGAMT(I) = MIN(CFAC*HEAT(I)/WSCALE(I),GAMCRT)
	651	HGAMQ(I) = MIN(CFAC*EVAP(I)/WSCALE(I),GAMCRQ)
	652	VPERT = HGAMT(I) + FVTHE1(I)HGAMQ(I)
	653	VPERT = MIN(VPERT,GAMCRT)
	654	THERMAL(I) = THERMAL(I) + MAX(VPERT,RZERO)
	655	HGAMT(I) = MAX(HGAMT(I),RZERO)
	656	HGAMQ(I) = MAX(HGAMQ(I),RZERO)
	657	ELSE
	658	PBLFLG(I) = .FALSE.
	659	ENDIF
	660	ENDDO
	661	!
	662	DO I = 1,IM
	663	IF(PBLFLG(I)) THEN
	664	KPBL(I) = 1
	665	HPBL(I) = ZI(I,2)
	666	ENDIF
	667	ENDDO
	668	!
	669	! ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL
	670	!
	671	DO I = 1, IM
	672	IF(PBLFLG(I)) THEN
	673	STABLE(I) = .FALSE.
	674	RBUP(I) = RBSOIL(I)
	675	ENDIF
	676	ENDDO
	677	DO K = 2, KMPBL
	678	DO I = 1, IM
	679	IF(.NOT.STABLE(I).AND.PBLFLG(I)) THEN
	680	RBDN(I) = RBUP(I)
	681	! ZL(I,k) = ZL(I,K-1) - (T1(i,k)+T1(i,K-1))/2 *
	682	! & LOG(PRSL(I,K)/PRSL(I,K-1)) * ROG
	683	THEKV(I) = THETA(i,k)(1.+FVMAX(Q1(i,k,1),QMIN))
	684	SPDK2 = MAX(RCL(i)(U1(i,k)U1(i,k)+V1(i,k)*V1(i,k)),RONE)
	685	RBUP(I) = (THEKV(I)-THERMAL(I))(GZL(I,k)/THE1V(I))/SPDK2
	686	KPBL(I) = K
	687	STABLE(I) = RBUP(I).GT.RBCR
	688	ENDIF
	689	ENDDO
	690	ENDDO
	691	!
	692	DO I = 1,IM
	693	IF(PBLFLG(I)) THEN
	694	K = KPBL(I)
	695	IF(RBDN(I).GE.RBCR) THEN
	696	RBINT = 0.
	697	ELSEIF(RBUP(I).LE.RBCR) THEN
	698	RBINT = 1.
	699	ELSE
	700	RBINT = (RBCR-RBDN(I))/(RBUP(I)-RBDN(I))
	701	ENDIF
	702	HPBL(I) = ZL(I,K-1) + RBINT*(ZL(I,k)-ZL(I,K-1))
	703	IF(HPBL(I).LT.ZI(I,KPBL(I))) KPBL(I) = KPBL(I) - 1
	704	IF(KPBL(I).LE.1) PBLFLG(I) = .FALSE.
	705	ENDIF
	706	ENDDO
	707	!!
	708	!
	709	! COMPUTE DIFFUSION COEFFICIENTS BELOW PBL
	710	!
	711	DO K = 1, KMPBL
	712	DO I=1,IM
	713	IF(KPBL(I).GT.K) THEN
	714	PRINV = 1.0 / (PHIH(I)/PHIM(I)+CFACVK.1)
	715	PRINV = MIN(PRINV,PRMAX)
	716	PRINV = MAX(PRINV,PRMIN)
	717	! ZFAC = MAX((1.-(ZI(I,K+1)-ZL1(I))/ &
	718	! & (HPBL(I)-ZL1(I))), ZFMIN)
	719	ZFAC = MAX((1.-(ZI(I,K+1)-ZL(I,1))/ &
	720	& (HPBL(I)-ZL(I,1))), ZFMIN)
	721	DKU(i,k) = XKZO + WSCALE(I)VKZI(I,K+1) &
	722	& * ZFAC**PFAC
	723	DKT(i,k) = DKU(i,k)*PRINV
	724	DKO(i,k) = (DKU(i,k)-XKZO)*PRINV
	725	DKU(i,k) = MIN(DKU(i,k),DKMAX)
	726	DKU(i,k) = MAX(DKU(i,k),DKMIN)
	727	DKT(i,k) = MIN(DKT(i,k),DKMAX)
	728	DKT(i,k) = MAX(DKT(i,k),DKMIN)
	729	DKO(i,k) = MAX(RZERO, MIN(DKMAX, DKO(i,k)))
	730	ENDIF
	731	ENDDO
	732	ENDDO
	733	!
	734	! COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE)
	735	!
	736	DO K = 1, KM-1
	737	DO I=1,IM
	738	IF(K.GE.KPBL(I)) THEN
	739	! TI = 0.5*(T1(i,k)+T1(i,K+1))
	740	TI = 2.0 / (T1(i,k)+T1(i,K+1))
	741	! RDZ = RDZT(I,K)/TI
	742	RDZ = RDZT(I,K) * TI
	743	! RDZ = RDZT(I,K)
	744	DW2 = RCL(i)((U1(i,k)-U1(i,K+1))*2 &
	745	& + (V1(i,k)-V1(i,K+1))**2)
	746	SHR2 = MAX(DW2,DW2MIN)RDZ*2
	747	TVD = T1(i,k)(1.+FVMAX(Q1(i,k,1),QMIN))
	748	TVU = T1(i,K+1)(1.+FVMAX(Q1(i,K+1,1),QMIN))
	749	! BVF2 = G(GOCP+RDZ(TVU-TVD))/TI
	750	BVF2 = G(GOCP+RDZ(TVU-TVD)) * TI
	751	RI = MAX(BVF2/SHR2,RIMIN)
	752	ZK = VK*ZI(I,K+1)
	753	! RL2 = (ZKRLAM/(RLAM+ZK))*2
	754	! DK = RL2*SQRT(SHR2)
	755	RL2 = ZK*RLAM/(RLAM+ZK)
	756	DK = RL2RL2SQRT(SHR2)
	757	IF(RI.LT.0.) THEN ! UNSTABLE REGIME
	758	SRI = SQRT(-RI)
	759	DKU(i,k) = XKZO + DK(1+8.(-RI)/(1+1.746*SRI))
	760	! DKT(i,k) = XKZO + DK(1+8.(-RI)/(1+1.286*SRI))
	761	tem = DK(1+8.(-RI)/(1+1.286*SRI))
	762	DKT(i,k) = XKZO + tem
	763	DKO(i,k) = tem
	764	ELSE ! STABLE REGIME
	765	! DKT(i,k) = XKZO + DK/(1+5.RI)*2
	766	tem = DK/(1+5.RI)*2
	767	DKT(i,k) = XKZO + tem
	768	DKO(i,k) = tem
	769	PRNUM = 1.0 + 2.1*RI
	770	PRNUM = MIN(PRNUM,PRMAX)
	771	DKU(i,k) = (DKT(i,k)-XKZO)*PRNUM + XKZO
	772	ENDIF
	773	!
	774	DKU(i,k) = MIN(DKU(i,k),DKMAX)
	775	DKU(i,k) = MAX(DKU(i,k),DKMIN)
	776	DKT(i,k) = MIN(DKT(i,k),DKMAX)
	777	DKT(i,k) = MAX(DKT(i,k),DKMIN)
	778	DKO(i,k) = MAX(RZERO, MIN(DKMAX, DKO(i,k)))
	779	!
	780	!!! IF(I.EQ.LOND.AND.LAT.EQ.LATD) THEN
	781	!!! PRNUM = DKU(k)/DKT(k)
	782	!!! WRITE(IUN,610) K,PRNUM,DKT(k),DKU(k),RL2,RI,
	783	!!! 1 BVF2,SHR2
	784	!!! ENDIF
	785	!
	786	ENDIF
	787	ENDDO
	788	ENDDO
	789	!
	790	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE
	791	!
	792	DO I=1,IM
	793	AD(I,1) = 1.
	794	A1(I,1) = T1(i,1) + BETA(i) * HEAT(I)
	795	A2(I,1) = Q1(i,1,1) + BETA(i) * EVAP(I)
	796	! A1(I,1) = T1(i,1)-BETAT(I)*(THETA(i,1)-TSEA(I))
	797	! A2(I,1) = Q1(i,1,1)-BETAQ(I)*
	798	! & (MAX(Q1(i,1,1),QMIN)-MAX(QSS(I),QMIN))
	799	ENDDO
	800	!
	801	DO K = 1,KM-1
	802	DO I = 1,IM
	803	DTODSD = DT/DEL(I,K)
	804	DTODSU = DT/DEL(I,K+1)
	805	DSIG = PRSL(I,K)-PRSL(I,K+1)
	806	RDZ = RDZT(I,K)*2./(T1(i,k)+T1(i,K+1))
	807	! RDZ = RDZT(I,K)
	808	tem1 = DSIG * DKT(i,k) * RDZ
	809	IF(PBLFLG(I).AND.K.LT.KPBL(I)) THEN
	810	! DSDZT = DSIGDKT(i,k)RDZ*(GOCP-HGAMT(I)/HPBL(I))
	811	! DSDZQ = DSIGDKT(i,k)RDZ*(-HGAMQ(I)/HPBL(I))
	812	tem = 1.0 / HPBL(I)
	813	DSDZT = tem1 * (GOCP-HGAMT(I)*tem)
	814	DSDZQ = tem1 * (-HGAMQ(I)*tem)
	815	A2(I,k) = A2(I,k)+DTODSD*DSDZQ
	816	A2(I,k+1) = Q1(i,k+1,1)-DTODSU*DSDZQ
	817	ELSE
	818	! DSDZT = DSIGDKT(i,k)RDZ*(GOCP)
	819	DSDZT = tem1 * GOCP
	820	A2(I,k+1) = Q1(i,k+1,1)
	821	ENDIF
	822	! DSDZ2 = DSIGDKT(i,k)RDZ*RDZ
	823	DSDZ2 = tem1 * RDZ
	824	AU(I,k) = -DTODSD*DSDZ2
	825	AL(I,k) = -DTODSU*DSDZ2
	826	AD(I,k) = AD(I,k)-AU(I,k)
	827	AD(I,k+1) = 1.-AL(I,k)
	828	A1(I,k) = A1(I,k)+DTODSD*DSDZT
	829	A1(I,k+1) = T1(i,k+1)-DTODSU*DSDZT
	830	ENDDO
	831	ENDDO
	832	!
	833	! SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE
	834	!
	835	CALL TRIDIN(IM,KM,1,AL,AD,AU,A1,A2,AU,A1,A2)
	836	!
	837	! RECOVER TENDENCIES OF HEAT AND MOISTURE
	838	!
	839	DO K = 1,KM
	840	DO I = 1,IM
	841	TTEND = (A1(I,k)-T1(i,k))*RDT
	842	QTEND = (A2(I,k)-Q1(i,k,1))*RDT
	843	TAU(i,k) = TAU(i,k)+TTEND
	844	RTG(I,k,1) = RTG(i,k,1)+QTEND
	845	DTSFC(I) = DTSFC(I)+CONTDEL(I,K)TTEND
	846	DQSFC(I) = DQSFC(I)+CONQDEL(I,K)QTEND
	847	ENDDO
	848	ENDDO
	849	!
	850	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM
	851	!
	852	DO I=1,IM
	853	! AD(I,1) = 1.+BETAW(I)
	854	AD(I,1) = 1.0 + BETA(i) * STRESS(I) / SPD1(I)
	855	A1(I,1) = U1(i,1)
	856	A2(I,1) = V1(i,1)
	857	! AD(I,1) = 1.0
	858	! tem = 1.0 + BETA(I) * STRESS(I) / SPD1(I)
	859	! A1(I,1) = U1(i,1) * tem
	860	! A2(I,1) = V1(i,1) * tem
	861	ENDDO
	862	!
	863	DO K = 1,KM-1
	864	DO I=1,IM
	865	DTODSD = DT/DEL(I,K)
	866	DTODSU = DT/DEL(I,K+1)
	867	DSIG = PRSL(I,K)-PRSL(I,K+1)
	868	RDZ = RDZT(I,K)*2./(T1(i,k)+T1(i,k+1))
	869	! RDZ = RDZT(I,K)
	870	DSDZ2 = DSIGDKU(i,k)RDZ*RDZ
	871	AU(I,k) = -DTODSD*DSDZ2
	872	AL(I,k) = -DTODSU*DSDZ2
	873	AD(I,k) = AD(I,k)-AU(I,k)
	874	AD(I,k+1) = 1.-AL(I,k)
	875	A1(I,k+1) = U1(i,k+1)
	876	A2(I,k+1) = V1(i,k+1)
	877	ENDDO
	878	ENDDO
	879	!
	880	! SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM
	881	!
	882	CALL TRIDI2(IM,KM,AL,AD,AU,A1,A2,AU,A1,A2)
	883	!
	884	! RECOVER TENDENCIES OF MOMENTUM
	885	!
	886	DO K = 1,KM
	887	DO I = 1,IM
	888	CONWRC = CONW*SQRT(RCL(i))
	889	UTEND = (A1(I,k)-U1(i,k))*RDT
	890	VTEND = (A2(I,k)-V1(i,k))*RDT
	891	DU(i,k) = DU(i,k)+UTEND
	892	DV(i,k) = DV(i,k)+VTEND
	893	DUSFC(I) = DUSFC(I)+CONWRCDEL(I,K)UTEND
	894	DVSFC(I) = DVSFC(I)+CONWRCDEL(I,K)VTEND
	895	ENDDO
	896	ENDDO
	897	!!
	898	!
	899	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR TRACERS
	900	!
	901	if (ntrac .ge. 2) then
	902	DO I=1,IM
	903	AD(I,1) = 1.
	904	ENDDO
	905	do k = 2, ntrac
	906	is = (k-2) * km
	907	do i = 1, im
	908	AT(I,1+is) = Q1(i,1,k)
	909	enddo
	910	enddo
	911	!
	912	DO K = 1,KM-1
	913	DO I = 1,IM
	914	DTODSD = DT/DEL(I,K)
	915	DTODSU = DT/DEL(I,K+1)
	916	DSIG = PRSL(I,K)-PRSL(I,K+1)
	917	RDZ = RDZT(I,K)*2./(T1(i,k)+T1(i,K+1))
	918	tem1 = DSIG * DKT(i,k) * RDZ
	919	DSDZ2 = tem1 * RDZ
	920	AU(I,k) = -DTODSD*DSDZ2
	921	AL(I,k) = -DTODSU*DSDZ2
	922	AD(I,k) = AD(I,k)-AU(I,k)
	923	AD(I,k+1) = 1.-AL(I,k)
	924	ENDDO
	925	ENDDO
	926	do kk = 2, ntrac
	927	is = (kk-2) * km
	928	do k = 1, km - 1
	929	do i = 1, im
	930	AT(I,k+1+is) = Q1(i,k+1,kk)
	931	enddo
	932	enddo
	933	enddo
	934	!
	935	! SOLVE TRIDIAGONAL PROBLEM FOR TRACERS
	936	!
	937	CALL TRIDIT(IM,KM,ntrac-1,AL,AD,AU,AT,AU,AT)
	938	!
	939	! RECOVER TENDENCIES OF TRACERS
	940	!
	941	do kk = 2, ntrac
	942	is = (kk-2) * km
	943	do k = 1, km
	944	do i = 1, im
	945	QTEND = (AT(I,K+is)-Q1(i,K,kk))*RDT
	946	RTG(i,K,kk) = RTG(i,K,kk) + QTEND
	947	enddo
	948	enddo
	949	enddo
	950	endif
	951	!!
	952	RETURN
	953	END SUBROUTINE MONINP
	954	!FPP$ NOCONCUR R
	955	!-----------------------------------------------------------------------
	956	SUBROUTINE TRIDI2(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
	957	!sela %INCLUDE DBTRIDI2;
	958	!
	959	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	960	implicit none
	961	integer k,n,l,i
	962	real(kind=kind_phys) fk
	963	!
	964	real(kind=kind_phys) CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), &
	965	& AU(L,N-1),A1(L,N),A2(L,N)
	966	!-----------------------------------------------------------------------
	967	DO I=1,L
	968	FK = 1./CM(I,1)
	969	AU(I,1) = FK*CU(I,1)
	970	A1(I,1) = FK*R1(I,1)
	971	A2(I,1) = FK*R2(I,1)
	972	ENDDO
	973	DO K=2,N-1
	974	DO I=1,L
	975	FK = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
	976	AU(I,K) = FK*CU(I,K)
	977	A1(I,K) = FK(R1(I,K)-CL(I,K)A1(I,K-1))
	978	A2(I,K) = FK(R2(I,K)-CL(I,K)A2(I,K-1))
	979	ENDDO
	980	ENDDO
	981	DO I=1,L
	982	FK = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
	983	A1(I,N) = FK(R1(I,N)-CL(I,N)A1(I,N-1))
	984	A2(I,N) = FK(R2(I,N)-CL(I,N)A2(I,N-1))
	985	ENDDO
	986	DO K=N-1,1,-1
	987	DO I=1,L
	988	A1(I,K) = A1(I,K)-AU(I,K)*A1(I,K+1)
	989	A2(I,K) = A2(I,K)-AU(I,K)*A2(I,K+1)
	990	ENDDO
	991	ENDDO
	992	!-----------------------------------------------------------------------
	993	RETURN
	994	END SUBROUTINE TRIDI2
	995	!FPP$ NOCONCUR R
	996	!-----------------------------------------------------------------------
	997	SUBROUTINE TRIDIN(L,N,nt,CL,CM,CU,R1,R2,AU,A1,A2)
	998	!sela %INCLUDE DBTRIDI2;
	999	!
	1000	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	1001	implicit none
	1002	integer is,k,kk,n,nt,l,i
	1003	real(kind=kind_phys) fk(L)
	1004	!
	1005	real(kind=kind_phys) CL(L,2:N), CM(L,N), CU(L,N-1), &
	1006	& R1(L,N), R2(L,N*nt), &
	1007	& AU(L,N-1), A1(L,N), A2(L,N*nt), &
	1008	& FKK(L,2:N-1)
	1009	!-----------------------------------------------------------------------
	1010	DO I=1,L
	1011	FK(I) = 1./CM(I,1)
	1012	AU(I,1) = FK(I)*CU(I,1)
	1013	A1(I,1) = FK(I)*R1(I,1)
	1014	ENDDO
	1015	do k = 1, nt
	1016	is = (k-1) * n
	1017	do i = 1, l
	1018	a2(i,1+is) = fk(I) * r2(i,1+is)
	1019	enddo
	1020	enddo
	1021	DO K=2,N-1
	1022	DO I=1,L
	1023	FKK(I,K) = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
	1024	AU(I,K) = FKK(I,K)*CU(I,K)
	1025	A1(I,K) = FKK(I,K)(R1(I,K)-CL(I,K)A1(I,K-1))
	1026	ENDDO
	1027	ENDDO
	1028	do kk = 1, nt
	1029	is = (kk-1) * n
	1030	DO K=2,N-1
	1031	DO I=1,L
	1032	A2(I,K+is) = FKK(I,K)(R2(I,K+is)-CL(I,K)A2(I,K+is-1))
	1033	ENDDO
	1034	ENDDO
	1035	ENDDO
	1036	DO I=1,L
	1037	FK(I) = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
	1038	A1(I,N) = FK(I)(R1(I,N)-CL(I,N)A1(I,N-1))
	1039	ENDDO
	1040	do k = 1, nt
	1041	is = (k-1) * n
	1042	do i = 1, l
	1043	A2(I,N+is) = FK(I)(R2(I,N+is)-CL(I,N)A2(I,N+is-1))
	1044	enddo
	1045	enddo
	1046	DO K=N-1,1,-1
	1047	DO I=1,L
	1048	A1(I,K) = A1(I,K) - AU(I,K)*A1(I,K+1)
	1049	ENDDO
	1050	ENDDO
	1051	do kk = 1, nt
	1052	is = (kk-1) * n
	1053	DO K=n-1,1,-1
	1054	DO I=1,L
	1055	A2(I,K+is) = A2(I,K+is) - AU(I,K)*A2(I,K+is+1)
	1056	ENDDO
	1057	ENDDO
	1058	ENDDO
	1059	!-----------------------------------------------------------------------
	1060	RETURN
	1061	END SUBROUTINE TRIDIN
	1062	SUBROUTINE TRIDIT(L,N,nt,CL,CM,CU,RT,AU,AT)
	1063	!sela %INCLUDE DBTRIDI2;
	1064	!
	1065	USE MODULE_GFS_MACHINE , ONLY : kind_phys
	1066	implicit none
	1067	integer is,k,kk,n,nt,l,i
	1068	real(kind=kind_phys) fk(L)
	1069	!
	1070	real(kind=kind_phys) CL(L,2:N), CM(L,N), CU(L,N-1), &
	1071	& RT(L,N*nt), &
	1072	& AU(L,N-1), AT(L,N*nt), &
	1073	& FKK(L,2:N-1)
	1074	!-----------------------------------------------------------------------
	1075	DO I=1,L
	1076	FK(I) = 1./CM(I,1)
	1077	AU(I,1) = FK(I)*CU(I,1)
	1078	ENDDO
	1079	do k = 1, nt
	1080	is = (k-1) * n
	1081	do i = 1, l
	1082	at(i,1+is) = fk(I) * rt(i,1+is)
	1083	enddo
	1084	enddo
	1085	DO K=2,N-1
	1086	DO I=1,L
	1087	FKK(I,K) = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
	1088	AU(I,K) = FKK(I,K)*CU(I,K)
	1089	ENDDO
	1090	ENDDO
	1091	do kk = 1, nt
	1092	is = (kk-1) * n
	1093	DO K=2,N-1
	1094	DO I=1,L
	1095	AT(I,K+is) = FKK(I,K)(RT(I,K+is)-CL(I,K)AT(I,K+is-1))
	1096	ENDDO
	1097	ENDDO
	1098	ENDDO
	1099	DO I=1,L
	1100	FK(I) = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
	1101	ENDDO
	1102	do k = 1, nt
	1103	is = (k-1) * n
	1104	do i = 1, l
	1105	AT(I,N+is) = FK(I)(RT(I,N+is)-CL(I,N)AT(I,N+is-1))
	1106	enddo
	1107	enddo
	1108	do kk = 1, nt
	1109	is = (kk-1) * n
	1110	DO K=n-1,1,-1
	1111	DO I=1,L
	1112	AT(I,K+is) = AT(I,K+is) - AU(I,K)*AT(I,K+is+1)
	1113	ENDDO
	1114	ENDDO
	1115	ENDDO
	1116	!-----------------------------------------------------------------------
	1117	RETURN
	1118	END SUBROUTINE TRIDIT
	1119
	1120	!-----------------------------------------------------------------------
	1121
	1122	END MODULE module_bl_gfs

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

Download in other formats:

Original Format