Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_slab.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: 33.5 KB

Rev	Line
[2759]	1	!WRF:MODEL_LAYER:PHYSICS
	2	!
	3	MODULE module_sf_slab
	4
	5	!---SPECIFY CONSTANTS AND LAYERS FOR SOIL MODEL
	6	!---SOIL DIFFUSION CONSTANT SET (M^2/S)
	7
	8	REAL, PARAMETER :: DIFSL=5.e-7
	9
	10	!---FACTOR TO MAKE SOIL STEP MORE CONSERVATIVE
	11
	12	REAL , PARAMETER :: SOILFAC=1.25
	13
	14	CONTAINS
	15
	16	!----------------------------------------------------------------
	17	SUBROUTINE SLAB(T3D,QV3D,P3D,FLHC,FLQC, &
	18	PSFC,XLAND,TMN,HFX,QFX,LH,TSK,QSFC,CHKLOWQ, &
	19	GSW,GLW,CAPG,THC,SNOWC,EMISS,MAVAIL, &
	20	DELTSM,ROVCP,XLV,DTMIN,IFSNOW, &
	21	SVP1,SVP2,SVP3,SVPT0,EP2, &
	22	KARMAN,EOMEG,STBOLT, &
	23	TSLB,ZS,DZS,num_soil_layers,radiation, &
	24	P1000mb, &
	25	ids,ide, jds,jde, kds,kde, &
	26	ims,ime, jms,jme, kms,kme, &
	27	its,ite, jts,jte, kts,kte, &
	28	tml,t0ml,hml,h0ml,huml,hvml,ust,u_phy,v_phy, &
	29	f,g,omlcall,oml_gamma )
	30	!----------------------------------------------------------------
	31	IMPLICIT NONE
	32	!----------------------------------------------------------------
	33	!
	34	! SUBROUTINE SLAB CALCULATES THE GROUND TEMPERATURE TENDENCY
	35	! ACCORDING TO THE RESIDUAL OF THE SURFACE ENERGY BUDGET
	36	! (BLACKADAR, 1978B).
	37	!
	38	! CHANGES:
	39	! FOR SOIL SUB-TIMESTEPS UPDATE SURFACE HFX AND QFX AS TG
	40	! CHANGES TO PREVENT POSSIBLE INSTABILITY FOR LONG MODEL
	41	! STEPS (DT > ~200 SEC).
	42	!
	43	! PUT SNOW COVER CHECK ON SOIL SUB-TIMESTEPS
	44	!
	45	! MAKE UPPER LIMIT ON SOIL SUB-STEP LENGTH MORE CONSERVATIVE
	46	!
	47	!----------------------------------------------------------------
	48	!-- T3D temperature (K)
	49	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
	50	!-- P3D 3D pressure (Pa)
	51	!-- FLHC exchange coefficient for heat (m/s)
	52	!-- FLQC exchange coefficient for moisture (m/s)
	53	!-- PSFC surface pressure (Pa)
	54	!-- XLAND land mask (1 for land, 2 for water)
	55	!-- TMN soil temperature at lower boundary (K)
	56	!-- HFX upward heat flux at the surface (W/m^2)
	57	!-- QFX upward moisture flux at the surface (kg/m^2/s)
	58	!-- LH latent heat flux at the surface (W/m^2)
	59	!-- TSK surface temperature (K)
	60	!-- GSW downward short wave flux at ground surface (W/m^2)
	61	!-- GLW downward long wave flux at ground surface (W/m^2)
	62	!-- CAPG heat capacity for soil (J/K/m^3)
	63	!-- THC thermal inertia (Cal/cm/K/s^0.5)
	64	!-- SNOWC flag indicating snow coverage (1 for snow cover)
	65	!-- EMISS surface emissivity (between 0 and 1)
	66	!-- DELTSM time step (second)
	67	!-- ROVCP R/CP
	68	!-- XLV latent heat of melting (J/kg)
	69	!-- DTMIN time step (minute)
	70	!-- IFSNOW ifsnow=1 for snow-cover effects
	71	!-- SVP1 constant for saturation vapor pressure (kPa)
	72	!-- SVP2 constant for saturation vapor pressure (dimensionless)
	73	!-- SVP3 constant for saturation vapor pressure (K)
	74	!-- SVPT0 constant for saturation vapor pressure (K)
	75	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
	76	!-- EP2 constant for specific humidity calculation
	77	! (R_d/R_v) (dimensionless)
	78	!-- KARMAN Von Karman constant
	79	!-- EOMEG angular velocity of earth's rotation (rad/s)
	80	!-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
	81	!-- TSLB soil temperature in 5-layer model
	82	!-- ZS depths of centers of soil layers
	83	!-- DZS thicknesses of soil layers
	84	!-- TML ocean mixed layer temperature (K)
	85	!-- T0ML ocean mixed layer temperature (K) at initial time
	86	!-- HML ocean mixed layer depth (m)
	87	!-- H0ML ocean mixed layer depth (m) at initial time
	88	!-- HUML ocean mixed layer u component of wind
	89	!-- HVML ocean mixed layer v component of wind
	90	!-- OML_GAMMA deep water lapse rate (K m-1)
	91	!-- OMLCALL whether to call oml model
	92	!-- num_soil_layers the number of soil layers
	93	!-- ids start index for i in domain
	94	!-- ide end index for i in domain
	95	!-- jds start index for j in domain
	96	!-- jde end index for j in domain
	97	!-- kds start index for k in domain
	98	!-- kde end index for k in domain
	99	!-- ims start index for i in memory
	100	!-- ime end index for i in memory
	101	!-- jms start index for j in memory
	102	!-- jme end index for j in memory
	103	!-- kms start index for k in memory
	104	!-- kme end index for k in memory
	105	!-- its start index for i in tile
	106	!-- ite end index for i in tile
	107	!-- jts start index for j in tile
	108	!-- jte end index for j in tile
	109	!-- kts start index for k in tile
	110	!-- kte end index for k in tile
	111	!----------------------------------------------------------------
	112	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	113	ims,ime, jms,jme, kms,kme, &
	114	its,ite, jts,jte, kts,kte
	115
	116	INTEGER, INTENT(IN) :: num_soil_layers
	117	LOGICAL, INTENT(IN) :: radiation
	118
	119	INTEGER, INTENT(IN ) :: IFSNOW
	120
	121	!
	122	REAL, INTENT(IN ) :: DTMIN,XLV,ROVCP,DELTSM
	123
	124	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
	125	REAL, INTENT(IN ) :: EP2,KARMAN,EOMEG,STBOLT
	126	REAL, INTENT(IN ) :: P1000mb
	127
	128	REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
	129	INTENT(INOUT) :: TSLB
	130
	131	REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::ZS,DZS
	132
	133	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	134	INTENT(IN ) :: QV3D, &
	135	P3D, &
	136	T3D
	137	!
	138	REAL, DIMENSION( ims:ime, jms:jme ) , &
	139	INTENT(IN ) :: SNOWC, &
	140	XLAND, &
	141	EMISS, &
	142	MAVAIL, &
	143	TMN, &
	144	GSW, &
	145	GLW, &
	146	THC
	147
	148	!CHKLOWQ is declared as memory size
	149	!
	150	REAL, DIMENSION( ims:ime, jms:jme ) , &
	151	INTENT(INOUT) :: HFX, &
	152	QFX, &
	153	LH, &
	154	CAPG, &
	155	TSK, &
	156	QSFC, &
	157	CHKLOWQ
	158
	159	REAL, DIMENSION( ims:ime, jms:jme ) , &
	160	INTENT(IN ) :: PSFC
	161	!
	162	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
	163	FLHC, &
	164	FLQC
	165
	166	! Ocean Mixed Layer Vars
	167	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) :: &
	168	TML,T0ML,HML,H0ML,HUML,HVML
	169	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, INTENT(IN ) :: &
	170	U_PHY,V_PHY
	171	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN ) :: &
	172	UST, F
	173	REAL, OPTIONAL, INTENT(IN ) :: G
	174	REAL, OPTIONAL, INTENT(IN ) :: OML_GAMMA
	175	INTEGER, OPTIONAL, INTENT(IN ) :: OMLCALL
	176
	177	! LOCAL VARS
	178
	179	REAL, DIMENSION( its:ite ) :: QV1D, &
	180	P1D, &
	181	T1D
	182	INTEGER :: I,J
	183
	184	DO J=jts,jte
	185
	186	DO i=its,ite
	187	T1D(i) =T3D(i,1,j)
	188	QV1D(i)=QV3D(i,1,j)
	189	P1D(i) =P3D(i,1,j)
	190	ENDDO
	191
	192	! the indices to the PSFC argument in the following call look
	193	! wrong; however, it is correct to call with its (and not ims)
	194	! because of the way PSFC is defined in SLAB1D. Whether that
	195	! is a good idea or not, this commenter cannot comment. JM
	196
	197	CALL SLAB1D(J,T1D,QV1D,P1D,FLHC(ims,j),FLQC(ims,j), &
	198	PSFC(its,j),XLAND(ims,j),TMN(ims,j),HFX(ims,j), &
	199	QFX(ims,j),TSK(ims,j),QSFC(ims,j),CHKLOWQ(ims,j), &
	200	LH(ims,j),GSW(ims,j),GLW(ims,j), &
	201	CAPG(ims,j),THC(ims,j),SNOWC(ims,j),EMISS(ims,j), &
	202	MAVAIL(ims,j),DELTSM,ROVCP,XLV,DTMIN,IFSNOW, &
	203	SVP1,SVP2,SVP3,SVPT0,EP2,KARMAN,EOMEG,STBOLT, &
	204	TSLB(ims,1,j),ZS,DZS,num_soil_layers,radiation, &
	205	P1000mb, &
	206	ids,ide, jds,jde, kds,kde, &
	207	ims,ime, jms,jme, kms,kme, &
	208	its,ite, jts,jte, kts,kte )
	209
	210	IF (OMLCALL .EQ. 1) THEN
	211	! CALL wrf_debug( 1, 'Call OCEANML' )
	212	IF (PRESENT(tml) .AND. PRESENT(t0ml) &
	213	.AND. .TRUE. ) THEN
	214	DO i=its,ite
	215	IF(XLAND(I,J).GT.1.5)THEN
	216	CALL OCEANML(I,J,TML(i,j),T0ML(i,j),HML(i,j),H0ML(i,j),&
	217	HUML(i,j),HVML(i,j),TSK(i,j),HFX(i,j), &
	218	LH(i,j),GSW(i,j),GLW(i,j), &
	219	U_PHY(i,kts,j),V_PHY(i,kts,j),UST(i,j),F(i,j), &
	220	EMISS(i,j),STBOLT,G,DELTSM,OML_GAMMA, &
	221	ids,ide, jds,jde, kds,kde, &
	222	ims,ime, jms,jme, kms,kme, &
	223	its,ite, jts,jte, kts,kte )
	224	ENDIF
	225	ENDDO
	226	ELSE
	227	CALL wrf_error_fatal('Lacking arguments for OCEANML in slab')
	228	ENDIF
	229	ENDIF
	230
	231	ENDDO
	232
	233	END SUBROUTINE SLAB
	234
	235	!----------------------------------------------------------------
	236	SUBROUTINE SLAB1D(J,T1D,QV1D,P1D,FLHC,FLQC, &
	237	PSFCPA,XLAND,TMN,HFX,QFX,TSK,QSFC,CHKLOWQ, &
	238	LH,GSW,GLW,CAPG,THC,SNOWC,EMISS,MAVAIL, &
	239	DELTSM,ROVCP,XLV,DTMIN,IFSNOW, &
	240	SVP1,SVP2,SVP3,SVPT0,EP2, &
	241	KARMAN,EOMEG,STBOLT, &
	242	TSLB2D,ZS,DZS,num_soil_layers,radiation, &
	243	P1000mb, &
	244	ids,ide, jds,jde, kds,kde, &
	245	ims,ime, jms,jme, kms,kme, &
	246	its,ite, jts,jte, kts,kte )
	247	!----------------------------------------------------------------
	248	IMPLICIT NONE
	249	!----------------------------------------------------------------
	250	!
	251	! SUBROUTINE SLAB CALCULATES THE GROUND TEMPERATURE TENDENCY
	252	! ACCORDING TO THE RESIDUAL OF THE SURFACE ENERGY BUDGET
	253	! (BLACKADAR, 1978B).
	254	!
	255	! CHANGES:
	256	! FOR SOIL SUB-TIMESTEPS UPDATE SURFACE HFX AND QFX AS TG
	257	! CHANGES TO PREVENT POSSIBLE INSTABILITY FOR LONG MODEL
	258	! STEPS (DT > ~200 SEC).
	259	!
	260	! PUT SNOW COVER CHECK ON SOIL SUB-TIMESTEPS
	261	!
	262	! MAKE UPPER LIMIT ON SOIL SUB-STEP LENGTH MORE CONSERVATIVE
	263	!
	264	!----------------------------------------------------------------
	265
	266	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	267	ims,ime, jms,jme, kms,kme, &
	268	its,ite, jts,jte, kts,kte,J
	269
	270	INTEGER , INTENT(IN) :: num_soil_layers
	271	LOGICAL, INTENT(IN ) :: radiation
	272
	273	INTEGER, INTENT(IN ) :: IFSNOW
	274	!
	275	REAL, INTENT(IN ) :: DTMIN,XLV,ROVCP,DELTSM
	276
	277	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
	278	REAL, INTENT(IN ) :: EP2,KARMAN,EOMEG,STBOLT
	279	REAL, INTENT(IN ) :: P1000mb
	280
	281	REAL, DIMENSION( ims:ime , 1:num_soil_layers ), &
	282	INTENT(INOUT) :: TSLB2D
	283
	284	REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::ZS,DZS
	285
	286	!
	287	REAL, DIMENSION( ims:ime ) , &
	288	INTENT(INOUT) :: HFX, &
	289	QFX, &
	290	LH, &
	291	CAPG, &
	292	TSK, &
	293	QSFC, &
	294	CHKLOWQ
	295	!
	296	REAL, DIMENSION( ims:ime ) , &
	297	INTENT(IN ) :: SNOWC, &
	298	XLAND, &
	299	EMISS, &
	300	MAVAIL, &
	301	TMN, &
	302	GSW, &
	303	GLW, &
	304	THC
	305	!
	306	REAL, DIMENSION( its:ite ) , &
	307	INTENT(IN ) :: QV1D, &
	308	P1D, &
	309	T1D
	310	!
	311	REAL, DIMENSION( its:ite ) , &
	312	INTENT(IN ) :: PSFCPA
	313
	314	!
	315	REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: &
	316	FLHC, &
	317	FLQC
	318	! LOCAL VARS
	319
	320	REAL, DIMENSION( its:ite ) :: PSFC
	321
	322	REAL, DIMENSION( its:ite ) :: &
	323	THX, &
	324	QX, &
	325	SCR3
	326
	327	REAL, DIMENSION( its:ite ) :: DTHGDT, &
	328	TG0, &
	329	THTMN, &
	330	XLD1, &
	331	TSCVN, &
	332	OLTG, &
	333	UPFLUX, &
	334	HM, &
	335	RNET, &
	336	XINET, &
	337	QS, &
	338	DTSDT
	339	!
	340	REAL, DIMENSION( its:ite, num_soil_layers ) :: FLUX
	341	!
	342	INTEGER :: I,K,NSOIL,ITSOIL,L,NK,RADSWTCH
	343	REAL :: PS,PS1,XLDCOL,TSKX,RNSOIL,RHOG1,RHOG2,RHOG3,LAMDAG
	344	REAL :: THG,ESG,QSG,HFXT,QFXT,CS,CSW,LAMG(4),THCON,PL
	345
	346	!----------------------------------------------------------------------
	347	!-----DETERMINE IF ANY POINTS IN COLUMN ARE LAND (RATHER THAN OCEAN)
	348	! POINTS. IF NOT, SKIP DOWN TO THE PRINT STATEMENTS SINCE OCEAN
	349	! SURFACE TEMPERATURES ARE NOT ALLOWED TO CHANGE.
	350	!
	351	! from sfcrad
	352	!----------------------------------------------------------------------
	353	DATA CSW/4.183E6/
	354	DATA LAMG/1.407E-8, -1.455E-5, 6.290E-3, 0.16857/
	355
	356	DO i=its,ite
	357	! in cmb
	358	PSFC(I)=PSFCPA(I)/1000.
	359	ENDDO
	360
	361
	362	DO I=its,ite
	363	! PL cmb
	364	PL=P1D(I)/1000.
	365	SCR3(I)=T1D(I)
	366	! THCON=(100./PL)**ROVCP
	367	THCON=(P1000mb0.001/PL)*ROVCP
	368	THX(I)=SCR3(I)*THCON
	369	QX(I)=0.
	370	ENDDO
	371
	372	! IF(IDRY.EQ.1) GOTO 81
	373	DO I=its,ite
	374	QX(I)=QV1D(I)
	375	ENDDO
	376	81 CONTINUE
	377
	378	!
	379	!-----THE SLAB THERMAL CAPACITY CAPG(I) ARE DEPENDENT ON:
	380	! THC(I) - SOIL THERMAL INERTIAL, ONLY.
	381	!
	382	DO I=its,ite
	383	CAPG(I)=3.298E6*THC(I)
	384	IF(num_soil_layers .gt. 1)THEN
	385
	386	! CAPG REPRESENTS SOIL HEAT CAPACITY (J/K/M^3) WHEN DIFSL=5.E-7 (M^2/S)
	387	! TO GIVE A CORRECT THERMAL INERTIA (=CAPG*DIFSL^0.5)
	388
	389	CAPG(I)=5.9114E7*THC(I)
	390	ENDIF
	391	ENDDO
	392	!
	393	XLDCOL=2.0
	394	DO 10 I=its,ite
	395	XLDCOL=AMIN1(XLDCOL,XLAND(I))
	396	10 CONTINUE
	397	!
	398	IF(XLDCOL.GT.1.5)GOTO 90
	399	!
	400	!
	401	!-----CONVERT SLAB TEMPERATURE TO POTENTIAL TEMPERATURE AND
	402	! SET XLD1(I) = 0. FOR OCEAN POINTS:
	403	!
	404	!
	405	DO 20 I=its,ite
	406	IF((XLAND(I)-1.5).GE.0)THEN
	407	XLD1(I)=0.
	408	ELSE
	409	XLD1(I)=1.
	410	ENDIF
	411	20 CONTINUE
	412	!
	413	!-----CONVERT 'TSK(THETAG)' TO 'TG' FOR 'IUP' CALCULATION ....
	414	! IF WE ARE USING THE BLACKADAR MULTI-LEVEL (HIGH-RESOLUTION)
	415	! PBL MODEL
	416	!
	417	DO 50 I=its,ite
	418	IF(XLD1(I).LT.0.5)GOTO 50
	419
	420	! PS cmb
	421	PS=PSFC(I)
	422
	423	! TSK is Temperature at gound sfc
	424	! TG0(I)=TSK(I)(PS0.01)**ROVCP
	425	TG0(I)=TSK(I)
	426	50 CONTINUE
	427	!
	428	!-----COMPUTE THE SURFACE ENERGY BUDGET:
	429	!
	430	! IF(ISOIL.EQ.1)NSOIL=1
	431	IF(num_soil_layers .gt. 1)NSOIL=1
	432
	433
	434	IF (radiation) then
	435	RADSWTCH=1
	436	ELSE
	437	RADSWTCH=0
	438	ENDIF
	439
	440	DO 70 I=its,ite
	441	IF(XLD1(I).LT.0.5)GOTO 70
	442	! OLTG(I)=TSK(I)(100./PSFC(I))*ROVCP
	443	OLTG(I)=TSK(I)(P1000mb0.001/PSFC(I))**ROVCP
	444	UPFLUX(I)=RADSWTCHSTBOLTTG0(I)**4
	445	XINET(I)=EMISS(I)*(GLW(I)-UPFLUX(I))
	446	RNET(I)=GSW(I)+XINET(I)
	447	HM(I)=1.18EOMEG(TG0(I)-TMN(I))
	448	! MOISTURE FLUX CALCULATED HERE (OVERWRITES SFC LAYER VALUE FOR LAND)
	449	ESG=SVP1EXP(SVP2(TG0(I)-SVPT0)/(TG0(I)-SVP3))
	450	QSG=EP2*ESG/(PSFC(I)-ESG)
	451	THG=TSK(I)(100./PSFC(I))*ROVCP
	452	HFX(I)=FLHC(I)*(THG-THX(I))
	453	QFX(I)=FLQC(I)*(QSG-QX(I))
	454	LH(I)=QFX(I)*XLV
	455	QS(I)=HFX(I)+QFX(I)*XLV
	456	! IF(ISOIL.EQ.0)THEN
	457	IF(num_soil_layers .EQ. 1)THEN
	458	DTHGDT(I)=(RNET(I)-QS(I))/CAPG(I)-HM(I)
	459	ELSE
	460	DTHGDT(I)=0.
	461	ENDIF
	462	70 CONTINUE
	463	! IF(ISOIL.EQ.1)THEN
	464	IF(num_soil_layers .gt. 1)THEN
	465	NSOIL=1+IFIX(SOILFAC4DIFSL/DZS(1)*DELTSM/DZS(1))
	466	RNSOIL=1./FLOAT(NSOIL)
	467	!
	468	! SOIL SUB-TIMESTEP
	469	!
	470	DO ITSOIL=1,NSOIL
	471	DO I=its,ite
	472	DO L=1,num_soil_layers-1
	473	IF(XLD1(I).LT.0.5)GOTO 75
	474	IF(L.EQ.1.AND.ITSOIL.GT.1)THEN
	475	! PS1=(PSFC(I)0.01)*ROVCP
	476	PS1=(PSFCPA(I)/P1000mb)**ROVCP
	477
	478	! for rk scheme A and B are the same
	479	PS=PSFC(I)
	480	THG=TSLB2D(I,1)/PS1
	481	ESG=SVP1EXP(SVP2(TSLB2D(I,1)-SVPT0)/(TSLB2D(I,1) &
	482	-SVP3))
	483	QSG=EP2*ESG/(PS-ESG)
	484	! UPDATE FLUXES FOR NEW GROUND TEMPERATURE
	485	HFXT=FLHC(I)*(THG-THX(I))
	486	QFXT=FLQC(I)*(QSG-QX(I))
	487	QS(I)=HFXT+QFXT*XLV
	488	! SUM HFX AND QFX OVER SOIL TIMESTEPS
	489	HFX(I)=HFX(I)+HFXT
	490	QFX(I)=QFX(I)+QFXT
	491	ENDIF
	492	FLUX(I,1)=RNET(I)-QS(I)
	493	FLUX(I,L+1)=-DIFSLCAPG(I)(TSLB2D(I,L+1)-TSLB2D(I,L))/( &
	494	ZS(L+1)-ZS(L))
	495	DTSDT(I)=-(FLUX(I,L+1)-FLUX(I,L))/(DZS(L)*CAPG(I))
	496	TSLB2D(I,L)=TSLB2D(I,L)+DTSDT(I)DELTSMRNSOIL
	497	IF(IFSNOW.EQ.1.AND.L.EQ.1)THEN
	498	IF((SNOWC(I).GT.0..AND.TSLB2D(I,1).GT.273.16))THEN
	499	TSLB2D(I,1)=273.16
	500	ENDIF
	501	ENDIF
	502	IF(L.EQ.1)DTHGDT(I)=DTHGDT(I)+RNSOIL*DTSDT(I)
	503	IF(ITSOIL.EQ.NSOIL.AND.L.EQ.1)THEN
	504	! AVERAGE HFX AND QFX OVER SOIL TIMESTEPS FOR OUTPUT TO PBL
	505	HFX(I)=HFX(I)*RNSOIL
	506	QFX(I)=QFX(I)*RNSOIL
	507	LH(I)=QFX(I)*XLV
	508	ENDIF
	509	75 CONTINUE
	510	ENDDO
	511	ENDDO
	512	ENDDO
	513	ENDIF
	514	!
	515	DO 80 I=its,ite
	516	IF(XLD1(I).LT.0.5) GOTO 80
	517	TSKX=TG0(I)+DELTSM*DTHGDT(I)
	518
	519	! TSK is temperature
	520	! TSK(I)=TSKX(100./PS1)*ROVCP
	521	TSK(I)=TSKX
	522	80 CONTINUE
	523
	524	!
	525	!-----MODIFY THE THE GROUND TEMPERATURE IF THE SNOW COVER EFFECTS ARE
	526	! CONSIDERED: LIMIT THE GROUND TEMPERATURE UNDER 0 C.
	527	!
	528	IF(IFSNOW.EQ.0)GOTO 90
	529	DO 85 I=its,ite
	530	IF(XLD1(I).LT.0.5)GOTO 85
	531	! PS1=(PSFC(I)0.01)*ROVCP
	532	! TSCVN(I)=TSK(I)*PS1
	533	TSCVN(I)=TSK(I)
	534	IF((SNOWC(I).GT.0..AND.TSCVN(I).GT.273.16))THEN
	535	TSCVN(I)=273.16
	536	ELSE
	537	TSCVN(I)=TSCVN(I)
	538	ENDIF
	539	! TSK(I)=TSCVN(I)/PS1
	540	TSK(I)=TSCVN(I)
	541	85 CONTINUE
	542	!
	543	90 CONTINUE
	544	DO I=its,ite
	545	! QSFC and CHKLOWQ needed by Eta PBL
	546	QSFC(I)=QX(I)+QFX(I)/FLQC(I)
	547	CHKLOWQ(I)=MAVAIL(I)
	548	ENDDO
	549	!
	550	140 CONTINUE
	551
	552	END SUBROUTINE SLAB1D
	553
	554	!================================================================
	555	SUBROUTINE slabinit(TSK,TMN, &
	556	TSLB,ZS,DZS,num_soil_layers, &
	557	allowed_to_read, start_of_simulation, &
	558	ids,ide, jds,jde, kds,kde, &
	559	ims,ime, jms,jme, kms,kme, &
	560	its,ite, jts,jte, kts,kte, &
	561	oml_hml0, omlcall, &
	562	tml,t0ml,hml,h0ml,huml,hvml )
	563	!----------------------------------------------------------------
	564	IMPLICIT NONE
	565	!----------------------------------------------------------------
	566	LOGICAL , INTENT(IN) :: allowed_to_read
	567	LOGICAL , INTENT(IN) :: start_of_simulation
	568	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	569	ims,ime, jms,jme, kms,kme, &
	570	its,ite, jts,jte, kts,kte
	571
	572	INTEGER, INTENT(IN ) :: num_soil_layers
	573	!
	574	REAL, DIMENSION( ims:ime , 1:num_soil_layers , jms:jme ), INTENT(INOUT) :: TSLB
	575
	576	REAL, DIMENSION(1:num_soil_layers), INTENT(IN) :: ZS,DZS
	577
	578	REAL, DIMENSION( ims:ime, jms:jme ) , &
	579	INTENT(IN) :: TSK, &
	580	TMN
	581	REAL, DIMENSION( ims:ime, jms:jme ) , OPTIONAL , &
	582	INTENT(INOUT) :: TML, T0ML, HML, H0ML, HUML, HVML
	583	REAL , OPTIONAL, INTENT(IN ) :: oml_hml0
	584	INTEGER, OPTIONAL, INTENT(IN ) :: omlcall
	585
	586	! LOCAR VAR
	587
	588	INTEGER :: L,J,I,itf,jtf
	589	CHARACTER*1024 message
	590
	591	!----------------------------------------------------------------
	592
	593	itf=min0(ite,ide-1)
	594	jtf=min0(jte,jde-1)
	595
	596	IF ( PRESENT(omlcall) .AND. omlcall .EQ. 1) THEN
	597	WRITE(message,*)'Initializing OML with HML0 = ', oml_hml0
	598	CALL wrf_debug (1, TRIM(message))
	599
	600	IF(start_of_simulation .AND. &
	601	PRESENT(tml) .AND. PRESENT(t0ml) ) THEN
	602	DO J=jts,jtf
	603	DO I=its,itf
	604	TML(I,J)=TSK(I,J)
	605	T0ML(I,J)=TSK(I,J)
	606	! MAY HAVE INPUT OF HML BUT FOR NOW SET HERE
	607	HML(I,J)=oml_hml0
	608	H0ML(I,J)=HML(I,J)
	609	HUML(I,J)=0.
	610	HVML(I,J)=0.
	611	ENDDO
	612	ENDDO
	613	ENDIF
	614	ENDIF
	615
	616	END SUBROUTINE slabinit
	617
	618	!-------------------------------------------------------------------
	619
	620	SUBROUTINE OCEANML(I,J,TML,T0ML,H,H0,HUML, &
	621	HVML,TSK,HFX, &
	622	LH,GSW,GLW, &
	623	UAIR,VAIR,UST,F,EMISS,STBOLT,G,DT,OML_GAMMA, &
	624	ids,ide, jds,jde, kds,kde, &
	625	ims,ime, jms,jme, kms,kme, &
	626	its,ite, jts,jte, kts,kte )
	627
	628	!----------------------------------------------------------------
	629	IMPLICIT NONE
	630	!----------------------------------------------------------------
	631	INTEGER, INTENT(IN ) :: I, J
	632	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	633	ims,ime, jms,jme, kms,kme, &
	634	its,ite, jts,jte, kts,kte
	635
	636	REAL, INTENT(INOUT) :: TML, H, H0, HUML, HVML, TSK
	637
	638	REAL, INTENT(IN ) :: T0ML, HFX, LH, GSW, GLW, &
	639	UAIR, VAIR, UST, F, EMISS
	640
	641	REAL, INTENT(IN) :: STBOLT, G, DT, OML_GAMMA
	642
	643	! Local
	644	REAL :: rhoair, rhowater, Gam, alp, BV2, A1, A2, B2, u, v, wspd, &
	645	hu1, hv1, hu2, hv2, taux, tauy, tauxair, tauyair, q, hold, &
	646	hsqrd, thp, cwater
	647	CHARACTER(LEN=120) :: time_series
	648
	649	hu1=huml
	650	hv1=hvml
	651	rhoair=1.
	652	rhowater=1000.
	653	cwater=4200.
	654	! Deep ocean lapse rate (K/m) - from Rich
	655	Gam=oml_gamma
	656	! if(i.eq.1 .and. j.eq.1 .or. i.eq.105.and.j.eq.105) print *, 'gamma = ', gam
	657	! Gam=0.14
	658	! Gam=5.6/40.
	659	! if(i.eq.1 .and. j.eq.1 ) print *, 'gamma = ', gam
	660	! Gam=5./100.
	661	! Thermal expansion coeff (/K)
	662	! alp=.0002
	663	! temp dependence (/K)
	664	alp=max((tml-273.15)*1.e-5, 1.e-6)
	665	BV2=alpgGam
	666	thp=t0ml-Gam*(h-h0)
	667	A1=(tml-thp)h - 0.5Gamhh
	668	if(h.ne.0.)then
	669	u=hu1/h
	670	v=hv1/h
	671	else
	672	u=0.
	673	v=0.
	674	endif
	675
	676	! time step
	677
	678	! q=(-hfx-lh+gsw+glw-stboltemisstmltmltmltml)/(rhowatercwater)
	679	! wspd=max(sqrt(uairuair+vairvair),0.1)
	680	wspd=sqrt(uairuair+vairvair)
	681	if (wspd .lt. 1.e-10 ) then
	682	! print *, 'i,j,wspd are ', i,j,wspd
	683	wspd = 1.e-10
	684	endif
	685	tauxair=ustustuair/wspd
	686	taux=rhoair/rhowater*tauxair
	687	tauyair=ustustvair/wspd
	688	tauy=rhoair/rhowater*tauyair
	689	! note: forward-backward coriolis force for effective time-centering
	690	hu2=hu1+dt( fhv1 + taux)
	691	hv2=hv1+dt(-fhu2 + tauy)
	692	! A2=A1+q*dt
	693	A2=A1
	694
	695	huml=hu2
	696	hvml=hv2
	697
	698	hold=h
	699	B2=hu2hu2+hv2hv2
	700	hsqrd=-A2/Gam + sqrt(A2A2/(GamGam) + 2.*B2/BV2)
	701	h=sqrt(max(hsqrd,0.0))
	702	! limit to positive h change
	703	if(h.lt.hold)h=hold
	704
	705	if(h.ne.0.)then
	706	tml=max(t0ml - Gam(h-h0) + 0.5Gam*h + A2/h, 273.15)
	707	u=hu2/h
	708	v=hv2/h
	709	else
	710	tml=t0ml
	711	u=0.
	712	v=0.
	713	endif
	714	tsk=tml
	715	! if(h.gt.100.)print *,i,j,h,tml,' h,tml'
	716
	717	! ww: output point data
	718	! if( (i.eq.190 .and. j.eq.115) .or. (i.eq.170 .and. j.eq.125) ) then
	719	! write(jtime,fmt='("TS ",f10.0)') float(itimestep)
	720	! CALL wrf_message ( TRIM(jtime) )
	721	! write(time_series,fmt='("OML",2I4,2F9.5,2F8.2,2E15.5,F8.3)') &
	722	! i,j,u,v,tml,h,taux,tauy,a2
	723	! CALL wrf_message ( TRIM(time_series) )
	724	! end if
	725
	726	END SUBROUTINE OCEANML
	727	!-------------------------------------------------------------------
	728
	729	END MODULE module_sf_slab

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

Download in other formats:

Original Format