Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_sfclay.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
File size: 45.1 KB

Rev	Line
[2759]	1	!WRF:MODEL_LAYER:PHYSICS
	2	!
	3	MODULE module_sf_sfclay
	4
	5	REAL , PARAMETER :: VCONVC=1.
	6	REAL , PARAMETER :: CZO=0.0185
	7	REAL , PARAMETER :: OZO=1.59E-5
	8
	9	REAL, DIMENSION(0:1000 ),SAVE :: PSIMTB,PSIHTB
	10
	11	CONTAINS
	12
	13	!-------------------------------------------------------------------
	14	SUBROUTINE SFCLAY(U3D,V3D,T3D,QV3D,P3D,dz8w, &
	15	CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
	16	ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
	17	XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
	18	U10,V10,TH2,T2,Q2, &
	19	GZ1OZ0,WSPD,BR,ISFFLX,DX, &
	20	SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
	21	KARMAN,EOMEG,STBOLT, &
	22	P1000mb, &
	23	ids,ide, jds,jde, kds,kde, &
	24	ims,ime, jms,jme, kms,kme, &
	25	its,ite, jts,jte, kts,kte, &
	26	uratx,vratx,tratx, &
	27	ustm,ck,cka,cd,cda,isftcflx )
	28	!-------------------------------------------------------------------
	29	IMPLICIT NONE
	30	!-------------------------------------------------------------------
	31	!-- U3D 3D u-velocity interpolated to theta points (m/s)
	32	!-- V3D 3D v-velocity interpolated to theta points (m/s)
	33	!-- T3D temperature (K)
	34	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
	35	!-- P3D 3D pressure (Pa)
	36	!-- dz8w dz between full levels (m)
	37	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
	38	!-- G acceleration due to gravity (m/s^2)
	39	!-- ROVCP R/CP
	40	!-- R gas constant for dry air (J/kg/K)
	41	!-- XLV latent heat of vaporization for water (J/kg)
	42	!-- PSFC surface pressure (Pa)
	43	!-- ZNT roughness length (m)
	44	!-- UST u* in similarity theory (m/s)
	45	!-- USTM u* in similarity theory (m/s) without vconv correction
	46	! used to couple with TKE scheme
	47	!-- PBLH PBL height from previous time (m)
	48	!-- MAVAIL surface moisture availability (between 0 and 1)
	49	!-- ZOL z/L height over Monin-Obukhov length
	50	!-- MOL T* (similarity theory) (K)
	51	!-- REGIME flag indicating PBL regime (stable, unstable, etc.)
	52	!-- PSIM similarity stability function for momentum
	53	!-- PSIH similarity stability function for heat
	54	!-- XLAND land mask (1 for land, 2 for water)
	55	!-- HFX upward heat flux at the surface (W/m^2)
	56	!-- QFX upward moisture flux at the surface (kg/m^2/s)
	57	!-- LH net upward latent heat flux at surface (W/m^2)
	58	!-- TSK surface temperature (K)
	59	!-- FLHC exchange coefficient for heat (W/m^2/K)
	60	!-- FLQC exchange coefficient for moisture (kg/m^2/s)
	61	!-- CHS heat/moisture exchange coefficient for LSM (m/s)
	62	!-- QGH lowest-level saturated mixing ratio
	63	!-- QSFC ground saturated mixing ratio
	64	!-- uratx ratio of surface U to U10
	65	!-- vratx ratio of surface V to V10
	66	!-- tratx ratio of surface T to TH2
	67	!-- U10 diagnostic 10m u wind
	68	!-- V10 diagnostic 10m v wind
	69	!-- TH2 diagnostic 2m theta (K)
	70	!-- T2 diagnostic 2m temperature (K)
	71	!-- Q2 diagnostic 2m mixing ratio (kg/kg)
	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	!-- ISFFLX isfflx=1 for surface heat and moisture fluxes
	76	!-- DX horizontal grid size (m)
	77	!-- SVP1 constant for saturation vapor pressure (kPa)
	78	!-- SVP2 constant for saturation vapor pressure (dimensionless)
	79	!-- SVP3 constant for saturation vapor pressure (K)
	80	!-- SVPT0 constant for saturation vapor pressure (K)
	81	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
	82	!-- EP2 constant for specific humidity calculation
	83	! (R_d/R_v) (dimensionless)
	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	!-- ck enthalpy exchange coeff at 10 meters
	88	!-- cd momentum exchange coeff at 10 meters
	89	!-- cka enthalpy exchange coeff at the lowest model level
	90	!-- cda momentum exchange coeff at the lowest model level
	91	!-- isftcflx =0, (Charnock and Carlson-Boland); =1, AHW Ck, Cd
	92	!-- ids start index for i in domain
	93	!-- ide end index for i in domain
	94	!-- jds start index for j in domain
	95	!-- jde end index for j in domain
	96	!-- kds start index for k in domain
	97	!-- kde end index for k in domain
	98	!-- ims start index for i in memory
	99	!-- ime end index for i in memory
	100	!-- jms start index for j in memory
	101	!-- jme end index for j in memory
	102	!-- kms start index for k in memory
	103	!-- kme end index for k in memory
	104	!-- its start index for i in tile
	105	!-- ite end index for i in tile
	106	!-- jts start index for j in tile
	107	!-- jte end index for j in tile
	108	!-- kts start index for k in tile
	109	!-- kte end index for k in tile
	110	!-------------------------------------------------------------------
	111	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	112	ims,ime, jms,jme, kms,kme, &
	113	its,ite, jts,jte, kts,kte
	114	!
	115	INTEGER, INTENT(IN ) :: ISFFLX
	116	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
	117	REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
	118	REAL, INTENT(IN ) :: P1000mb
	119	!
	120	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	121	INTENT(IN ) :: dz8w
	122
	123	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	124	INTENT(IN ) :: QV3D, &
	125	P3D, &
	126	T3D
	127
	128	REAL, DIMENSION( ims:ime, jms:jme ) , &
	129	INTENT(IN ) :: MAVAIL, &
	130	PBLH, &
	131	XLAND, &
	132	TSK
	133	REAL, DIMENSION( ims:ime, jms:jme ) , &
	134	INTENT(OUT ) :: U10, &
	135	V10, &
	136	TH2, &
	137	T2, &
	138	Q2, &
	139	QSFC
	140
	141	!
	142	REAL, DIMENSION( ims:ime, jms:jme ) , &
	143	INTENT(INOUT) :: REGIME, &
	144	HFX, &
	145	QFX, &
	146	LH, &
	147	MOL,RMOL
	148	!m the following 5 are change to memory size
	149	!
	150	REAL, DIMENSION( ims:ime, jms:jme ) , &
	151	INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
	152	PSIM,PSIH
	153
	154	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	155	INTENT(IN ) :: U3D, &
	156	V3D
	157
	158	REAL, DIMENSION( ims:ime, jms:jme ) , &
	159	INTENT(IN ) :: PSFC
	160
	161	REAL, DIMENSION( ims:ime, jms:jme ) , &
	162	INTENT(INOUT) :: ZNT, &
	163	ZOL, &
	164	UST, &
	165	CPM, &
	166	CHS2, &
	167	CQS2, &
	168	CHS
	169
	170	REAL, DIMENSION( ims:ime, jms:jme ) , &
	171	INTENT(INOUT) :: FLHC,FLQC
	172
	173	REAL, DIMENSION( ims:ime, jms:jme ) , &
	174	INTENT(INOUT) :: &
	175	QGH
	176
	177
	178
	179	REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
	180
	181	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , &
	182	INTENT(OUT) :: uratx,vratx,tratx
	183
	184	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , &
	185	INTENT(OUT) :: ck,cka,cd,cda,ustm
	186
	187	INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX
	188
	189	! LOCAL VARS
	190
	191	REAL, DIMENSION( its:ite ) :: U1D, &
	192	V1D, &
	193	QV1D, &
	194	P1D, &
	195	T1D
	196
	197	REAL, DIMENSION( its:ite ) :: dz8w1d
	198
	199	INTEGER :: I,J
	200
	201	DO J=jts,jte
	202	DO i=its,ite
	203	dz8w1d(I) = dz8w(i,1,j)
	204	ENDDO
	205
	206	DO i=its,ite
	207	U1D(i) =U3D(i,1,j)
	208	V1D(i) =V3D(i,1,j)
	209	QV1D(i)=QV3D(i,1,j)
	210	P1D(i) =P3D(i,1,j)
	211	T1D(i) =T3D(i,1,j)
	212	ENDDO
	213
	214	! Sending array starting locations of optional variables may cause
	215	! troubles, so we explicitly change the call.
	216
	217	CALL SFCLAY1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d, &
	218	CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
	219	CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j), &
	220	ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j), &
	221	MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j), &
	222	XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j), &
	223	U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j), &
	224	Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j), &
	225	QSFC(ims,j),LH(ims,j), &
	226	GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX, &
	227	SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT, &
	228	P1000mb, &
	229	ids,ide, jds,jde, kds,kde, &
	230	ims,ime, jms,jme, kms,kme, &
	231	its,ite, jts,jte, kts,kte &
	232	#if ( EM_CORE == 1 )
	233	,uratx(ims,j),vratx(ims,j),tratx(ims,j), &
	234	isftcflx, &
	235	USTM(ims,j),CK(ims,j),CKA(ims,j), &
	236	CD(ims,j),CDA(ims,j) &
	237	#endif
	238	)
	239	ENDDO
	240
	241
	242	END SUBROUTINE SFCLAY
	243
	244
	245	!-------------------------------------------------------------------
	246	SUBROUTINE SFCLAY1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d, &
	247	CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
	248	ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
	249	XLAND,HFX,QFX,TSK, &
	250	U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH, &
	251	QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, &
	252	SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
	253	KARMAN,EOMEG,STBOLT, &
	254	P1000mb, &
	255	ids,ide, jds,jde, kds,kde, &
	256	ims,ime, jms,jme, kms,kme, &
	257	its,ite, jts,jte, kts,kte, &
	258	uratx,vratx,tratx, &
	259	isftcflx, &
	260	ustm,ck,cka,cd,cda )
	261	!-------------------------------------------------------------------
	262	IMPLICIT NONE
	263	!-------------------------------------------------------------------
	264	REAL, PARAMETER :: XKA=2.4E-5
	265	REAL, PARAMETER :: PRT=1.
	266
	267	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	268	ims,ime, jms,jme, kms,kme, &
	269	its,ite, jts,jte, kts,kte, &
	270	J
	271	!
	272	INTEGER, INTENT(IN ) :: ISFFLX
	273	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
	274	REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
	275	REAL, INTENT(IN ) :: P1000mb
	276
	277	!
	278	REAL, DIMENSION( ims:ime ) , &
	279	INTENT(IN ) :: MAVAIL, &
	280	PBLH, &
	281	XLAND, &
	282	TSK
	283	!
	284	REAL, DIMENSION( ims:ime ) , &
	285	INTENT(IN ) :: PSFCPA
	286
	287	REAL, DIMENSION( ims:ime ) , &
	288	INTENT(INOUT) :: REGIME, &
	289	HFX, &
	290	QFX, &
	291	MOL,RMOL
	292	!m the following 5 are changed to memory size---
	293	!
	294	REAL, DIMENSION( ims:ime ) , &
	295	INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
	296	PSIM,PSIH
	297
	298	REAL, DIMENSION( ims:ime ) , &
	299	INTENT(INOUT) :: ZNT, &
	300	ZOL, &
	301	UST, &
	302	CPM, &
	303	CHS2, &
	304	CQS2, &
	305	CHS
	306
	307	REAL, DIMENSION( ims:ime ) , &
	308	INTENT(INOUT) :: FLHC,FLQC
	309
	310	REAL, DIMENSION( ims:ime ) , &
	311	INTENT(INOUT) :: &
	312	QGH
	313
	314	REAL, DIMENSION( ims:ime ) , &
	315	INTENT(OUT) :: U10,V10, &
	316	TH2,T2,Q2,QSFC,LH
	317
	318
	319	REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
	320
	321	! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
	322	REAL, DIMENSION( its:ite ), INTENT(IN ) :: dz8w1d
	323
	324	REAL, DIMENSION( its:ite ), INTENT(IN ) :: UX, &
	325	VX, &
	326	QV1D, &
	327	P1D, &
	328	T1D
	329
	330	REAL, OPTIONAL, DIMENSION( ims:ime ) , &
	331	INTENT(OUT) :: uratx,vratx,tratx
	332
	333	REAL, OPTIONAL, DIMENSION( ims:ime ) , &
	334	INTENT(OUT) :: ck,cka,cd,cda,ustm
	335
	336	INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX
	337
	338	! LOCAL VARS
	339
	340	REAL, DIMENSION( its:ite ) :: ZA, &
	341	THVX,ZQKL, &
	342	ZQKLP1, &
	343	THX,QX, &
	344	PSIH2, &
	345	PSIM2, &
	346	PSIH10, &
	347	PSIM10, &
	348	DENOMQ, &
	349	WSPDI, &
	350	GZ2OZ0, &
	351	GZ10OZ0
	352	!
	353	REAL, DIMENSION( its:ite ) :: &
	354	RHOX,GOVRTH, &
	355	TGDSA
	356	!
	357	REAL, DIMENSION( its:ite) :: SCR3,SCR4
	358	REAL, DIMENSION( its:ite ) :: THGB, PSFC
	359	!
	360	INTEGER :: KL
	361
	362	INTEGER :: N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
	363
	364	REAL :: PL,THCON,TVCON,E1
	365	REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
	366	REAL :: DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
	367	REAL :: FLUXC,VSGD,Z0Q
	368	!-------------------------------------------------------------------
	369	KL=kte
	370
	371	DO i=its,ite
	372	! PSFC cb
	373	PSFC(I)=PSFCPA(I)/1000.
	374	ENDDO
	375	!
	376	!----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
	377	!
	378	DO 5 I=its,ite
	379	TGDSA(I)=TSK(I)
	380	! PSFC cb
	381	! THGB(I)=TSK(I)(100./PSFC(I))*ROVCP
	382	THGB(I)=TSK(I)(P1000mb/PSFCPA(I))*ROVCP
	383	5 CONTINUE
	384	!
	385	!-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
	386	! T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
	387	!
	388	! * NOTE *
	389	! THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
	390	! SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
	391	! TENDENCIES.
	392	!
	393	10 CONTINUE
	394
	395	! DO 24 I=its,ite
	396	! UX(I)=U1D(I)
	397	! VX(I)=V1D(I)
	398	! 24 CONTINUE
	399
	400	26 CONTINUE
	401
	402	!.....SCR3(I,K) STORE TEMPERATURE,
	403	! SCR4(I,K) STORE VIRTUAL TEMPERATURE.
	404
	405	DO 30 I=its,ite
	406	! PL cb
	407	PL=P1D(I)/1000.
	408	SCR3(I)=T1D(I)
	409	! THCON=(100./PL)**ROVCP
	410	THCON=(P1000mb0.001/PL)*ROVCP
	411	THX(I)=SCR3(I)*THCON
	412	SCR4(I)=SCR3(I)
	413	THVX(I)=THX(I)
	414	QX(I)=0.
	415	30 CONTINUE
	416	!
	417	DO I=its,ite
	418	QGH(I)=0.
	419	FLHC(I)=0.
	420	FLQC(I)=0.
	421	CPM(I)=CP
	422	ENDDO
	423	!
	424	! IF(IDRY.EQ.1)GOTO 80
	425	DO 50 I=its,ite
	426	QX(I)=QV1D(I)
	427	TVCON=(1.+EP1*QX(I))
	428	THVX(I)=THX(I)*TVCON
	429	SCR4(I)=SCR3(I)*TVCON
	430	50 CONTINUE
	431	!
	432	DO 60 I=its,ite
	433	E1=SVP1EXP(SVP2(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
	434	QSFC(I)=EP2*E1/(PSFC(I)-E1)
	435	! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
	436	! Q2SAT = QGH IN LSM
	437	E1=SVP1EXP(SVP2(T1D(I)-SVPT0)/(T1D(I)-SVP3))
	438	PL=P1D(I)/1000.
	439	QGH(I)=EP2*E1/(PL-E1)
	440	CPM(I)=CP(1.+0.8QX(I))
	441	60 CONTINUE
	442	80 CONTINUE
	443
	444	!-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
	445	! LEVEL, AND THE LAYER THICKNESSES.
	446
	447	DO 90 I=its,ite
	448	ZQKLP1(I)=0.
	449	RHOX(I)=PSFC(I)1000./(RSCR4(I))
	450	90 CONTINUE
	451	!
	452	DO 110 I=its,ite
	453	ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
	454	110 CONTINUE
	455	!
	456	DO 120 I=its,ite
	457	ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))
	458	120 CONTINUE
	459	!
	460	DO 160 I=its,ite
	461	GOVRTH(I)=G/THX(I)
	462	160 CONTINUE
	463
	464	!-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
	465	! AKB(1976), EQ(12).
	466
	467	DO 260 I=its,ite
	468	GZ1OZ0(I)=ALOG(ZA(I)/ZNT(I))
	469	GZ2OZ0(I)=ALOG(2./ZNT(I))
	470	GZ10OZ0(I)=ALOG(10./ZNT(I))
	471	IF((XLAND(I)-1.5).GE.0)THEN
	472	ZL=ZNT(I)
	473	ELSE
	474	ZL=0.01
	475	ENDIF
	476	WSPD(I)=SQRT(UX(I)UX(I)+VX(I)VX(I))
	477
	478	TSKV=THGB(I)(1.+EP1QSFC(I)*MAVAIL(I))
	479	DTHVDZ=(THVX(I)-TSKV)
	480	! Convective velocity scale Vc and subgrid-scale velocity Vsg
	481	! following Beljaars (1995, QJRMS) and Mahrt and Sun (1995, MWR)
	482	! ... HONG Aug. 2001
	483	!
	484	! VCONV = 0.25sqrt(g/tskvpblh(i)*dthvm)
	485	! Use Beljaars over land, old MM5 (Wyngaard) formula over water
	486	if (xland(i).lt.1.5) then
	487	fluxc = max(hfx(i)/rhox(i)/cp &
	488	+ ep1tskvqfx(i)/rhox(i),0.)
	489	VCONV = vconvc(g/tgdsa(i)pblh(i)fluxc)*.33
	490	else
	491	IF(-DTHVDZ.GE.0)THEN
	492	DTHVM=-DTHVDZ
	493	ELSE
	494	DTHVM=0.
	495	ENDIF
	496	VCONV = 2.*SQRT(DTHVM)
	497	endif
	498	! Mahrt and Sun low-res correction
	499	VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
	500	WSPD(I)=SQRT(WSPD(I)WSPD(I)+VCONVVCONV+vsgd*vsgd)
	501	WSPD(I)=AMAX1(WSPD(I),0.1)
	502	BR(I)=GOVRTH(I)ZA(I)DTHVDZ/(WSPD(I)*WSPD(I))
	503	! IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2
	504	IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
	505	!jdf
	506	RMOL(I)=-GOVRTH(I)DTHVDZZA(I)*KARMAN
	507	!jdf
	508
	509	260 CONTINUE
	510
	511	!
	512	!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
	513	!
	514	!
	515	! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
	516	! AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
	517	!
	518	! CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
	519	!
	520	! 1. BR .GE. 0.2;
	521	! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
	522	!
	523	! 2. BR .LT. 0.2 .AND. BR .GT. 0.0;
	524	! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
	525	! (REGIME=2),
	526	!
	527	! 3. BR .EQ. 0.0
	528	! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
	529	!
	530	! 4. BR .LT. 0.0
	531	! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
	532	!
	533	!CCCCC
	534
	535	DO 320 I=its,ite
	536	!CCCCC
	537	!CC REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
	538	!CC IF(BR(I).LT.0..AND.HOL(I,J).GT.1.5)GOTO 310
	539	IF(BR(I).LT.0.)GOTO 310
	540	!
	541	!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
	542	!
	543	IF(BR(I).LT.0.2)GOTO 270
	544	REGIME(I)=1.
	545	PSIM(I)=-10.*GZ1OZ0(I)
	546	! LOWER LIMIT ON PSI IN STABLE CONDITIONS
	547	PSIM(I)=AMAX1(PSIM(I),-10.)
	548	PSIH(I)=PSIM(I)
	549	PSIM10(I)=10./ZA(I)*PSIM(I)
	550	PSIM10(I)=AMAX1(PSIM10(I),-10.)
	551	PSIH10(I)=PSIM10(I)
	552	PSIM2(I)=2./ZA(I)*PSIM(I)
	553	PSIM2(I)=AMAX1(PSIM2(I),-10.)
	554	PSIH2(I)=PSIM2(I)
	555
	556	! 1.0 over Monin-Obukhov length
	557	IF(UST(I).LT.0.01)THEN
	558	RMOL(I)=BR(I)*GZ1OZ0(I) !ZA/L
	559	ELSE
	560	RMOL(I)=KARMANGOVRTH(I)ZA(I)MOL(I)/(UST(I)UST(I)) !ZA/L
	561	ENDIF
	562	RMOL(I)=AMIN1(RMOL(I),9.999) ! ZA/L
	563	RMOL(I) = RMOL(I)/ZA(I) !1.0/L
	564
	565	GOTO 320
	566	!
	567	!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
	568	!
	569	270 IF(BR(I).EQ.0.0)GOTO 280
	570	REGIME(I)=2.
	571	PSIM(I)=-5.0BR(I)GZ1OZ0(I)/(1.1-5.0*BR(I))
	572	! LOWER LIMIT ON PSI IN STABLE CONDITIONS
	573	PSIM(I)=AMAX1(PSIM(I),-10.)
	574	!.....AKB(1976), EQ(16).
	575	PSIH(I)=PSIM(I)
	576	PSIM10(I)=10./ZA(I)*PSIM(I)
	577	PSIM10(I)=AMAX1(PSIM10(I),-10.)
	578	PSIH10(I)=PSIM10(I)
	579	PSIM2(I)=2./ZA(I)*PSIM(I)
	580	PSIM2(I)=AMAX1(PSIM2(I),-10.)
	581	PSIH2(I)=PSIM2(I)
	582
	583	! Linear form: PSIM = -0.5*ZA/L; e.g, see eqn 16 of
	584	! Blackadar, Modeling the nocturnal boundary layer, Preprints,
	585	! Third Symposium on Atmospheric Turbulence Diffusion and Air Quality,
	586	! Raleigh, NC, 1976
	587	ZOL(I) = BR(I)GZ1OZ0(I)/(1.00001-5.0BR(I))
	588
	589	if ( ZOL(I) .GT. 0.5 ) then ! linear form ok
	590	! Holtslag and de Bruin, J. App. Meteor 27, 689-704, 1988;
	591	! see also, Launiainen, Boundary-Layer Meteor 76,165-179, 1995
	592	! Eqn (8) of Launiainen, 1995
	593	ZOL(I) = ( 1.89GZ1OZ0(I) + 44.2 ) BR(I)*BR(I) &
	594	+ ( 1.18GZ1OZ0(I) - 1.37 ) BR(I)
	595	ZOL(I)=AMIN1(ZOL(I),9.999)
	596	end if
	597
	598	! 1.0 over Monin-Obukhov length
	599	RMOL(I)= ZOL(I)/ZA(I)
	600
	601	GOTO 320
	602	!
	603	!-----CLASS 3; FORCED CONVECTION:
	604	!
	605	280 REGIME(I)=3.
	606	PSIM(I)=0.0
	607	PSIH(I)=PSIM(I)
	608	PSIM10(I)=0.
	609	PSIH10(I)=PSIM10(I)
	610	PSIM2(I)=0.
	611	PSIH2(I)=PSIM2(I)
	612
	613
	614	IF(UST(I).LT.0.01)THEN
	615	ZOL(I)=BR(I)*GZ1OZ0(I)
	616	ELSE
	617	ZOL(I)=KARMANGOVRTH(I)ZA(I)MOL(I)/(UST(I)UST(I))
	618	ENDIF
	619
	620	RMOL(I) = ZOL(I)/ZA(I)
	621
	622	GOTO 320
	623	!
	624	!-----CLASS 4; FREE CONVECTION:
	625	!
	626	310 CONTINUE
	627	REGIME(I)=4.
	628	IF(UST(I).LT.0.01)THEN
	629	ZOL(I)=BR(I)*GZ1OZ0(I)
	630	ELSE
	631	ZOL(I)=KARMANGOVRTH(I)ZA(I)MOL(I)/(UST(I)UST(I))
	632	ENDIF
	633	ZOL10=10./ZA(I)*ZOL(I)
	634	ZOL2=2./ZA(I)*ZOL(I)
	635	ZOL(I)=AMIN1(ZOL(I),0.)
	636	ZOL(I)=AMAX1(ZOL(I),-9.9999)
	637	ZOL10=AMIN1(ZOL10,0.)
	638	ZOL10=AMAX1(ZOL10,-9.9999)
	639	ZOL2=AMIN1(ZOL2,0.)
	640	ZOL2=AMAX1(ZOL2,-9.9999)
	641	NZOL=INT(-ZOL(I)*100.)
	642	RZOL=-ZOL(I)*100.-NZOL
	643	NZOL10=INT(-ZOL10*100.)
	644	RZOL10=-ZOL10*100.-NZOL10
	645	NZOL2=INT(-ZOL2*100.)
	646	RZOL2=-ZOL2*100.-NZOL2
	647	PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
	648	PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
	649	PSIM10(I)=PSIMTB(NZOL10)+RZOL10*(PSIMTB(NZOL10+1)-PSIMTB(NZOL10))
	650	PSIH10(I)=PSIHTB(NZOL10)+RZOL10*(PSIHTB(NZOL10+1)-PSIHTB(NZOL10))
	651	PSIM2(I)=PSIMTB(NZOL2)+RZOL2*(PSIMTB(NZOL2+1)-PSIMTB(NZOL2))
	652	PSIH2(I)=PSIHTB(NZOL2)+RZOL2*(PSIHTB(NZOL2+1)-PSIHTB(NZOL2))
	653
	654	!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
	655	!--- THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
	656	! PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
	657	! PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
	658	PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
	659	PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
	660	PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
	661	PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
	662	! AHW: mods to compute ck, cd
	663	PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I))
	664
	665	RMOL(I) = ZOL(I)/ZA(I)
	666
	667	320 CONTINUE
	668	!
	669	!-----COMPUTE THE FRICTIONAL VELOCITY:
	670	! ZA(1982) EQS(2.60),(2.61).
	671	!
	672	DO 330 I=its,ite
	673	DTG=THX(I)-THGB(I)
	674	PSIX=GZ1OZ0(I)-PSIM(I)
	675	PSIX10=GZ10OZ0(I)-PSIM10(I)
	676	! LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
	677	! ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
	678	PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)
	679
	680	IF((XLAND(I)-1.5).GE.0)THEN
	681	ZL=ZNT(I)
	682	ELSE
	683	ZL=0.01
	684	ENDIF
	685	PSIQ=ALOG(KARMANUST(I)ZA(I)/XKA+ZA(I)/ZL)-PSIH(I)
	686	PSIT2=GZ2OZ0(I)-PSIH2(I)
	687	PSIQ2=ALOG(KARMANUST(I)2./XKA+2./ZL)-PSIH2(I)
	688	! AHW: mods to compute ck, cd
	689	PSIQ10=ALOG(KARMANUST(I)10./XKA+10./ZL)-PSIH10(I)
	690	IF ( PRESENT(ISFTCFLX) ) THEN
	691	IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN
	692	Z0Q = 1.e-4 + 1.e-3(MAX(0.,UST(I)-1.))*2
	693	! Z0Q = 0.621.5E-5/UST(I) + 1.E-3(MAX(0.,UST(I)-1.))**2
	694	PSIQ=ALOG(ZA(I)/Z0Q)-PSIH(I)
	695	PSIT=PSIQ
	696	PSIQ2=ALOG(2./Z0Q)-PSIH2(I)
	697	PSIQ10=ALOG(10./Z0Q)-PSIH10(I)
	698	PSIT2=PSIQ2
	699	ENDIF
	700	ENDIF
	701	IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN
	702	Ck(I)=(karman/psix10)*(karman/psiq10)
	703	Cd(I)=(karman/psix10)*(karman/psix10)
	704	Cka(I)=(karman/psix)*(karman/psiq)
	705	Cda(I)=(karman/psix)*(karman/psix)
	706	ENDIF
	707	! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
	708	UST(I)=0.5UST(I)+0.5KARMAN*WSPD(I)/PSIX
	709	! TKE coupling: compute ust without vconv for use in tke scheme
	710	WSPDI(I)=SQRT(UX(I)UX(I)+VX(I)VX(I))
	711	IF ( PRESENT(USTM) ) THEN
	712	USTM(I)=0.5USTM(I)+0.5KARMAN*WSPDI(I)/PSIX
	713	ENDIF
	714	U10(I)=UX(I)*PSIX10/PSIX
	715	V10(I)=VX(I)*PSIX10/PSIX
	716	TH2(I)=THGB(I)+DTG*PSIT2/PSIT
	717	Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ
	718	! T2(I) = TH2(I)(PSFC(I)/100.)*ROVCP
	719	T2(I) = TH2(I)(PSFCPA(I)/P1000mb)*ROVCP
	720	! LATER Q2 WILL BE OVERWRITTEN FOR LAND POINTS IN SURFCE
	721	! QA2(I,J) = Q2(I)
	722	! UA10(I,J) = U10(I)
	723	! VA10(I,J) = V10(I)
	724	! write(,1002)UST(I),KARMANWSPD(I),PSIX,KARMAN*WSPD(I)/PSIX
	725	!
	726	IF(PRESENT(uratx) .and. PRESENT(vratx) .and. PRESENT(tratx))THEN
	727	IF(ABS(U10(I)) .GT. 1.E-10) THEN
	728	uratx(I) = UX(I)/U10(I)
	729	ELSE
	730	uratx(I) = 1.2
	731	END IF
	732	IF(ABS(V10(I)) .GT. 1.E-10) THEN
	733	vratx(I) = VX(I)/V10(I)
	734	ELSE
	735	vratx(I) = 1.2
	736	END IF
	737	tratx(I) = THX(I)/TH2(I)
	738	ENDIF
	739
	740	IF((XLAND(I)-1.5).LT.0.)THEN
	741	UST(I)=AMAX1(UST(I),0.1)
	742	ENDIF
	743	MOL(I)=KARMAN*DTG/PSIT/PRT
	744	DENOMQ(I)=PSIQ
	745	330 CONTINUE
	746	!
	747	335 CONTINUE
	748
	749	!-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:
	750
	751	DO i=its,ite
	752	QFX(i)=0.
	753	HFX(i)=0.
	754	ENDDO
	755
	756	IF (ISFFLX.EQ.0) GOTO 410
	757
	758	!-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).
	759
	760	DO 360 I=its,ite
	761	IF((XLAND(I)-1.5).GE.0)THEN
	762	ZNT(I)=CZOUST(I)UST(I)/G+OZO
	763	! AHW: change roughness length, and hence the drag coefficients Ck and Cd
	764	IF ( PRESENT(ISFTCFLX) ) THEN
	765	IF ( ISFTCFLX.EQ.1 ) THEN
	766	ZNT(I)=10.exp(-9.UST(I)**(-.3333))
	767	ZNT(I)=MIN(ZNT(I),2.85e-3)
	768	ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
	769	! ZNT(I)=MAX(ZNT(I),1.27e-7)
	770	ENDIF
	771	ENDIF
	772	ZL = ZNT(I)
	773	ELSE
	774	ZL = 0.01
	775	ENDIF
	776	FLQC(I)=RHOX(I)MAVAIL(I)UST(I)*KARMAN/DENOMQ(I)
	777	! FLQC(I)=RHOX(I)MAVAIL(I)UST(I)*KARMAN/( &
	778	! ALOG(KARMANUST(I)ZA(I)/XKA+ZA(I)/ZL)-PSIH(I))
	779	DTTHX=ABS(THX(I)-THGB(I))
	780	IF(DTTHX.GT.1.E-5)THEN
	781	FLHC(I)=CPM(I)RHOX(I)UST(I)*MOL(I)/(THX(I)-THGB(I))
	782	! write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I
	783	1001 format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
	784	ELSE
	785	FLHC(I)=0.
	786	ENDIF
	787	360 CONTINUE
	788
	789	!
	790	!-----COMPUTE SURFACE MOIST FLUX:
	791	!
	792	! IF(IDRY.EQ.1)GOTO 390
	793	!
	794	DO 370 I=its,ite
	795	QFX(I)=FLQC(I)*(QSFC(I)-QX(I))
	796	QFX(I)=AMAX1(QFX(I),0.)
	797	LH(I)=XLV*QFX(I)
	798	370 CONTINUE
	799
	800	!-----COMPUTE SURFACE HEAT FLUX:
	801	!
	802	390 CONTINUE
	803	DO 400 I=its,ite
	804	IF(XLAND(I)-1.5.GT.0.)THEN
	805	HFX(I)=FLHC(I)*(THGB(I)-THX(I))
	806	ELSEIF(XLAND(I)-1.5.LT.0.)THEN
	807	HFX(I)=FLHC(I)*(THGB(I)-THX(I))
	808	HFX(I)=AMAX1(HFX(I),-250.)
	809	ENDIF
	810	400 CONTINUE
	811
	812	DO I=its,ite
	813	IF((XLAND(I)-1.5).GE.0)THEN
	814	ZL=ZNT(I)
	815	ELSE
	816	ZL=0.01
	817	ENDIF
	818	CHS(I)=UST(I)KARMAN/(ALOG(KARMANUST(I)*ZA(I) &
	819	/XKA+ZA(I)/ZL)-PSIH(I))
	820	! GZ2OZ0(I)=ALOG(2./ZNT(I))
	821	! PSIM2(I)=-10.*GZ2OZ0(I)
	822	! PSIM2(I)=AMAX1(PSIM2(I),-10.)
	823	! PSIH2(I)=PSIM2(I)
	824	CQS2(I)=UST(I)KARMAN/(ALOG(KARMANUST(I)*2.0 &
	825	/XKA+2.0/ZL)-PSIH2(I))
	826	CHS2(I)=UST(I)*KARMAN/(GZ2OZ0(I)-PSIH2(I))
	827	ENDDO
	828
	829	410 CONTINUE
	830	!jdf
	831	! DO I=its,ite
	832	! IF(UST(I).GE.0.1) THEN
	833	! RMOL(I)=RMOL(I)(-FLHC(I))/(UST(I)UST(I)*UST(I))
	834	! ELSE
	835	! RMOL(I)=RMOL(I)(-FLHC(I))/(0.10.1*0.1)
	836	! ENDIF
	837	! ENDDO
	838	!jdf
	839
	840	!
	841	END SUBROUTINE SFCLAY1D
	842
	843	!====================================================================
	844	SUBROUTINE sfclayinit( allowed_to_read )
	845
	846	LOGICAL , INTENT(IN) :: allowed_to_read
	847	INTEGER :: N
	848	REAL :: ZOLN,X,Y
	849
	850	DO N=0,1000
	851	ZOLN=-FLOAT(N)*0.01
	852	X=(1-16.ZOLN)*0.25
	853	PSIMTB(N)=2ALOG(0.5(1+X))+ALOG(0.5(1+XX))- &
	854	2.ATAN(X)+2.ATAN(1.)
	855	Y=(1-16ZOLN)*0.5
	856	PSIHTB(N)=2ALOG(0.5(1+Y))
	857	ENDDO
	858
	859	END SUBROUTINE sfclayinit
	860
	861	!-------------------------------------------------------------------
	862
	863	END MODULE module_sf_sfclay

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

Download in other formats:

Original Format