Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_initialize_real.F @ 2759

Last change on this file since 2759 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: 150.0 KB

Rev	Line
[2759]	1	!REAL:MODEL_LAYER:INITIALIZATION
	2
	3	! This MODULE holds the routines which are used to perform various initializations
	4	! for the individual domains, specifically for the Eulerian, mass-based coordinate.
	5
	6	!-----------------------------------------------------------------------
	7
	8	MODULE module_initialize_real
	9
	10	USE module_bc
	11	USE module_configure
	12	USE module_domain
	13	USE module_io_domain
	14	USE module_model_constants
	15	! USE module_si_io_nmm
	16	USE module_state_description
	17	USE module_timing
	18	USE module_soil_pre
	19	#ifdef DM_PARALLEL
	20	USE module_dm
	21	USE module_ext_internal
	22	#endif
	23
	24	INTEGER :: internal_time_loop
	25	INTEGER:: MPI_COMM_COMP,MYPE
	26
	27	CONTAINS
	28
	29	!-------------------------------------------------------------------
	30
	31	SUBROUTINE init_domain ( grid )
	32
	33	IMPLICIT NONE
	34
	35	! Input space and data. No gridded meteorological data has been stored, though.
	36
	37	! TYPE (domain), POINTER :: grid
	38	TYPE (domain) :: grid
	39
	40	! Local data.
	41
	42	INTEGER :: idum1, idum2
	43
	44	CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
	45
	46	CALL init_domain_nmm (grid &
	47	!
	48	#include <actual_args.inc>
	49	!
	50	)
	51
	52	END SUBROUTINE init_domain
	53
	54	!-------------------------------------------------------------------
	55	!---------------------------------------------------------------------
	56	SUBROUTINE init_domain_nmm ( grid &
	57	!
	58	# include <dummy_args.inc>
	59	!
	60	)
	61
	62	USE module_optional_input
	63	IMPLICIT NONE
	64
	65	! Input space and data. No gridded meteorological data has been stored, though.
	66
	67	! TYPE (domain), POINTER :: grid
	68	TYPE (domain) :: grid
	69
	70	# include <dummy_decl.inc>
	71
	72	TYPE (grid_config_rec_type) :: config_flags
	73
	74	! Local domain indices and counters.
	75
	76	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
	77
	78	INTEGER :: &
	79	ids, ide, jds, jde, kds, kde, &
	80	ims, ime, jms, jme, kms, kme, &
	81	its, ite, jts, jte, kts, kte, &
	82	ips, ipe, jps, jpe, kps, kpe, &
	83	i, j, k, NNXP, NNYP
	84
	85	! Local data
	86
	87	CHARACTER(LEN=19):: start_date
	88
	89	#ifdef DM_PARALLEL
	90
	91	LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
	92
	93	! INTEGER :: DOMDESC
	94	REAL,ALLOCATABLE :: SICE_G(:,:), SM_G(:,:)
	95	INTEGER, ALLOCATABLE:: IHE_G(:),IHW_G(:)
	96	INTEGER, ALLOCATABLE:: ITEMP(:,:)
	97	INTEGER :: my_e,my_n,my_s,my_w,my_ne,my_nw,my_se,my_sw,myi,myj,npe
	98	INTEGER :: istat,inpes,jnpes
	99	#endif
	100
	101
	102	CHARACTER (LEN=255) :: message
	103
	104	INTEGER :: error
	105	REAL :: p_surf, p_level
	106	REAL :: cof1, cof2
	107	REAL :: qvf , qvf1 , qvf2 , pd_surf
	108	REAL :: p00 , t00 , a
	109	REAL :: hold_znw, rmin,rmax
	110
	111	REAL :: p_top_requested , ptsgm
	112	INTEGER :: num_metgrid_levels, ICOUNT
	113	REAL , DIMENSION(max_eta) :: eta_levels
	114
	115	LOGICAL :: stretch_grid, dry_sounding, debug, log_flag_sst, hyb_coor
	116
	117	REAL, ALLOCATABLE,DIMENSION(:,:):: ADUM2D,SNOWC,HT,TG_ALT, &
	118	PDVP,PSFC_OUTV
	119
	120	REAL, ALLOCATABLE,DIMENSION(:,:,:):: P3D_OUT,P3DV_OUT,P3DV_IN, &
	121	QTMP,QTMP2
	122
	123	INTEGER, ALLOCATABLE, DIMENSION(:):: KHL2,KVL2,KHH2,KVH2, &
	124	KHLA,KHHA,KVLA,KVHA
	125
	126	! INTEGER, ALLOCATABLE, DIMENSION(:,:):: LU_INDEX
	127
	128	REAL, ALLOCATABLE, DIMENSION(:):: DXJ,WPDARJ,CPGFUJ,CURVJ, &
	129	FCPJ,FDIVJ,EMJ,EMTJ,FADJ, &
	130	HDACJ,DDMPUJ,DDMPVJ
	131
	132	REAL, ALLOCATABLE,DIMENSION(:),SAVE:: SG1,SG2,DSG1,DSG2, &
	133	SGML1,SGML2
	134
	135	!-- Carsel and Parrish [1988]
	136	REAL , DIMENSION(100) :: lqmi
	137	integer iicount
	138
	139	REAL:: TPH0D,TLM0D
	140	REAL:: TPH0,WB,SB,TDLM,TDPH
	141	REAL:: WBI,SBI,EBI,ANBI,STPH0,CTPH0
	142	REAL:: TSPH,DTAD,DTCF
	143	REAL:: ACDT,CDDAMP,DXP,FP
	144	REAL:: WBD,SBD
	145	REAL:: RSNOW,SNOFAC
	146	REAL, PARAMETER:: SALP=2.60
	147	REAL, PARAMETER:: SNUP=0.040
	148	REAL:: SMCSUM,STCSUM,SEAICESUM,FISX
	149	REAL:: cur_smc, aposs_smc
	150
	151	INTEGER,PARAMETER:: DOUBLE=SELECTED_REAL_KIND(15,300)
	152	REAL(KIND=DOUBLE):: TERM1,APH,TLM,TPH,DLM,DPH,STPH,CTPH
	153
	154	INTEGER:: KHH,KVH,JAM,JA, IHL, IHH, L
	155	INTEGER:: II,JJ,ISRCH,ISUM,ITER,Ilook,Jlook
	156
	157	REAL, PARAMETER:: DTR=0.01745329
	158	REAL, PARAMETER:: W_NMM=0.08
	159	REAL, PARAMETER:: COAC=1.6
	160	REAL, PARAMETER:: CODAMP=6.4
	161	REAL, PARAMETER:: TWOM=.00014584
	162	REAL, PARAMETER:: CP=1004.6
	163	REAL, PARAMETER:: DFC=1.0
	164	REAL, PARAMETER:: DDFC=8.0
	165	REAL, PARAMETER:: ROI=916.6
	166	REAL, PARAMETER:: R=287.04
	167	REAL, PARAMETER:: CI=2060.0
	168	REAL, PARAMETER:: ROS=1500.
	169	REAL, PARAMETER:: CS=1339.2
	170	REAL, PARAMETER:: DS=0.050
	171	REAL, PARAMETER:: AKS=.0000005
	172	REAL, PARAMETER:: DZG=2.85
	173	REAL, PARAMETER:: DI=.1000
	174	REAL, PARAMETER:: AKI=0.000001075
	175	REAL, PARAMETER:: DZI=2.0
	176	REAL, PARAMETER:: THL=210.
	177	REAL, PARAMETER:: PLQ=70000.
	178	REAL, PARAMETER:: ERAD=6371200.
	179	REAL, PARAMETER:: TG0=258.16
	180	REAL, PARAMETER:: TGA=30.0
	181
	182	if (ALLOCATED(ADUM2D)) DEALLOCATE(ADUM2D)
	183	if (ALLOCATED(TG_ALT)) DEALLOCATE(TG_ALT)
	184
	185	#define COPY_IN
	186	#include <scalar_derefs.inc>
	187	#ifdef DM_PARALLEL
	188	# include <data_calls.inc>
	189	#endif
	190
	191	SELECT CASE ( model_data_order )
	192	CASE ( DATA_ORDER_ZXY )
	193	kds = grid%sd31 ; kde = grid%ed31 ;
	194	ids = grid%sd32 ; ide = grid%ed32 ;
	195	jds = grid%sd33 ; jde = grid%ed33 ;
	196
	197	kms = grid%sm31 ; kme = grid%em31 ;
	198	ims = grid%sm32 ; ime = grid%em32 ;
	199	jms = grid%sm33 ; jme = grid%em33 ;
	200
	201	kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
	202	its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
	203	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
	204
	205	CASE ( DATA_ORDER_XYZ )
	206	ids = grid%sd31 ; ide = grid%ed31 ;
	207	jds = grid%sd32 ; jde = grid%ed32 ;
	208	kds = grid%sd33 ; kde = grid%ed33 ;
	209
	210	ims = grid%sm31 ; ime = grid%em31 ;
	211	jms = grid%sm32 ; jme = grid%em32 ;
	212	kms = grid%sm33 ; kme = grid%em33 ;
	213
	214	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
	215	jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
	216	kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
	217
	218	CASE ( DATA_ORDER_XZY )
	219	ids = grid%sd31 ; ide = grid%ed31 ;
	220	kds = grid%sd32 ; kde = grid%ed32 ;
	221	jds = grid%sd33 ; jde = grid%ed33 ;
	222
	223	ims = grid%sm31 ; ime = grid%em31 ;
	224	kms = grid%sm32 ; kme = grid%em32 ;
	225	jms = grid%sm33 ; jme = grid%em33 ;
	226
	227	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
	228	kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
	229	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
	230
	231	END SELECT
	232
	233	#ifdef DM_PARALLEL
	234	CALL WRF_GET_MYPROC(MYPE)
	235	CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
	236	call wrf_get_nprocx(inpes)
	237	call wrf_get_nprocy(jnpes)
	238	!
	239	allocate(itemp(inpes,jnpes),stat=istat)
	240	npe=0
	241	!
	242	do j=1,jnpes
	243	do i=1,inpes
	244	itemp(i,j)=npe
	245	if(npe==mype)then
	246	myi=i
	247	myj=j
	248	endif
	249	npe=npe+1
	250	enddo
	251	enddo
	252	!
	253	my_n=-1
	254	if(myj+1<=jnpes)my_n=itemp(myi,myj+1)
	255	!
	256	my_e=-1
	257	if(myi+1<=inpes)my_e=itemp(myi+1,myj)
	258	!
	259	my_s=-1
	260	if(myj-1>=1)my_s=itemp(myi,myj-1)
	261	!
	262	my_w=-1
	263	if(myi-1>=1)my_w=itemp(myi-1,myj)
	264	!
	265	my_ne=-1
	266	if((myi+1<=inpes).and.(myj+1<=jnpes)) &
	267	my_ne=itemp(myi+1,myj+1)
	268	!
	269	my_se=-1
	270	if((myi+1<=inpes).and.(myj-1>=1)) &
	271	my_se=itemp(myi+1,myj-1)
	272	!
	273	my_sw=-1
	274	if((myi-1>=1).and.(myj-1>=1)) &
	275	my_sw=itemp(myi-1,myj-1)
	276	!
	277	my_nw=-1
	278	if((myi-1>=1).and.(myj+1<=jnpes)) &
	279	my_nw=itemp(myi-1,myj+1)
	280	!
	281	! my_neb(1)=my_n
	282	! my_neb(2)=my_e
	283	! my_neb(3)=my_s
	284	! my_neb(4)=my_w
	285	! my_neb(5)=my_ne
	286	! my_neb(6)=my_se
	287	! my_neb(7)=my_sw
	288	! my_neb(8)=my_nw
	289	!
	290	deallocate(itemp)
	291	#endif
	292
	293	grid%DT=float(grid%TIME_STEP)
	294
	295	NNXP=min(ITE,IDE-1)
	296	NNYP=min(JTE,JDE-1)
	297
	298	write(message,*) 'IDE, JDE: ', IDE,JDE
	299	write(message,*) 'NNXP, NNYP: ', NNXP,NNYP
	300	CALL wrf_message(message)
	301
	302	JAM=6+2*(JDE-JDS-10)
	303
	304	if (internal_time_loop .eq. 1) then
	305	ALLOCATE(ADUM2D(grid%sm31:grid%em31,jms:jme))
	306	ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
	307	ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
	308	ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),&
	309	FADJ(JTS:NNYP))
	310	ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
	311	ALLOCATE(KHLA(JAM),KHHA(JAM))
	312	ALLOCATE(KVLA(JAM),KVHA(JAM))
	313	endif
	314
	315
	316	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
	317
	318	write(message,*) 'cen_lat: ', config_flags%cen_lat
	319	CALL wrf_debug(100,message)
	320	write(message,*) 'cen_lon: ', config_flags%cen_lon
	321	CALL wrf_debug(100,message)
	322	write(message,*) 'dx: ', config_flags%dx
	323	CALL wrf_debug(100,message)
	324	write(message,*) 'dy: ', config_flags%dy
	325	CALL wrf_debug(100,message)
	326	write(message,*) 'config_flags%start_year: ', config_flags%start_year
	327	CALL wrf_debug(100,message)
	328	write(message,*) 'config_flags%start_month: ', config_flags%start_month
	329	CALL wrf_debug(100,message)
	330	write(message,*) 'config_flags%start_day: ', config_flags%start_day
	331	CALL wrf_debug(100,message)
	332	write(message,*) 'config_flags%start_hour: ', config_flags%start_hour
	333	CALL wrf_debug(100,message)
	334
	335	write(start_date,435) config_flags%start_year, config_flags%start_month, &
	336	config_flags%start_day, config_flags%start_hour
	337	435 format(I4,'-',I2.2,'-',I2.2,'_',I2.2,':00:00')
	338
	339	dlmd=config_flags%dx
	340	dphd=config_flags%dy
	341	tph0d=config_flags%cen_lat
	342	tlm0d=config_flags%cen_lon
	343
	344	!==========================================================================
	345
	346	!!
	347
	348	! Check to see if the boundary conditions are set
	349	! properly in the namelist file.
	350	! This checks for sufficiency and redundancy.
	351
	352	CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
	353
	354	! Some sort of "this is the first time" initialization. Who knows.
	355
	356	grid%itimestep=0
	357
	358	! Pull in the info in the namelist to compare it to the input data.
	359
	360	grid%real_data_init_type = model_config_rec%real_data_init_type
	361	write(message,*) 'what is flag_metgrid: ', flag_metgrid
	362	CALL wrf_message(message)
	363
	364	IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
	365
	366	num_metgrid_levels = grid%num_metgrid_levels
	367
	368	IF (ght_gc(its,jts,num_metgrid_levels/2) .lt. ght_gc(its,jts,num_metgrid_levels/2+1)) THEN
	369
	370
	371	write(message,*) 'normal ground up file order'
	372	hyb_coor=.false.
	373	CALL wrf_message(message)
	374
	375	ELSE
	376
	377	hyb_coor=.true.
	378	write(message,*) 'reverse the order of coordinate'
	379	CALL wrf_message(message)
	380
	381	CALL reverse_vert_coord(ght_gc, 2, num_metgrid_levels &
	382	&, IDS,IDE,JDS,JDE,KDS,KDE &
	383	&, IMS,IME,JMS,JME,KMS,KME &
	384	&, ITS,ITE,JTS,JTE,KTS,KTE )
	385
	386	CALL reverse_vert_coord(p_gc, 2, num_metgrid_levels &
	387	&, IDS,IDE,JDS,JDE,KDS,KDE &
	388	&, IMS,IME,JMS,JME,KMS,KME &
	389	&, ITS,ITE,JTS,JTE,KTS,KTE )
	390
	391	CALL reverse_vert_coord(t_gc, 2, num_metgrid_levels &
	392	&, IDS,IDE,JDS,JDE,KDS,KDE &
	393	&, IMS,IME,JMS,JME,KMS,KME &
	394	&, ITS,ITE,JTS,JTE,KTS,KTE )
	395
	396	CALL reverse_vert_coord(u_gc, 2, num_metgrid_levels &
	397	&, IDS,IDE,JDS,JDE,KDS,KDE &
	398	&, IMS,IME,JMS,JME,KMS,KME &
	399	&, ITS,ITE,JTS,JTE,KTS,KTE )
	400
	401	CALL reverse_vert_coord(v_gc, 2, num_metgrid_levels &
	402	&, IDS,IDE,JDS,JDE,KDS,KDE &
	403	&, IMS,IME,JMS,JME,KMS,KME &
	404	&, ITS,ITE,JTS,JTE,KTS,KTE )
	405
	406	CALL reverse_vert_coord(rh_gc, 2, num_metgrid_levels &
	407	&, IDS,IDE,JDS,JDE,KDS,KDE &
	408	&, IMS,IME,JMS,JME,KMS,KME &
	409	&, ITS,ITE,JTS,JTE,KTS,KTE )
	410
	411	endif
	412
	413
	414	IF (hyb_coor) THEN
	415	! limit extreme deviations from source model topography
	416	! due to potential for nasty extrapolation/interpolation issues
	417	!
	418	write(message,*) 'min, max of ht_gc before adjust: ', minval(ht_gc), maxval(ht_gc)
	419	CALL wrf_debug(100,message)
	420	ICOUNT=0
	421	DO J=JTS,min(JTE,JDE-1)
	422	DO I=ITS,min(ITE,IDE-1)
	423	IF ((ht_gc(I,J) - ght_gc(I,J,2)) .LT. -150.) THEN
	424	ht_gc(I,J)=ght_gc(I,J,2)-150.
	425	IF (ICOUNT .LT. 20) THEN
	426	write(message,*) 'increasing NMM topo toward RUC ', I,J
	427	CALL wrf_debug(100,message)
	428	ICOUNT=ICOUNT+1
	429	ENDIF
	430	ELSEIF ((ht_gc(I,J) - ght_gc(I,J,2)) .GT. 150.) THEN
	431	ht_gc(I,J)=ght_gc(I,J,2)+150.
	432	IF (ICOUNT .LT. 20) THEN
	433	write(message,*) 'decreasing NMM topo toward RUC ', I,J
	434	CALL wrf_debug(100,message)
	435	ICOUNT=ICOUNT+1
	436	ENDIF
	437	ENDIF
	438	END DO
	439	END DO
	440
	441	write(message,*) 'min, max of ht_gc after correct: ', minval(ht_gc), maxval(ht_gc)
	442	CALL wrf_debug(100,message)
	443	ENDIF
	444
	445	CALL boundary_smooth(ht_gc,landmask, grid, 12 , 12 &
	446	&, IDS,IDE,JDS,JDE,KDS,KDE &
	447	&, IMS,IME,JMS,JME,KMS,KME &
	448	&, ITS,ITE,JTS,JTE,KTS,KTE )
	449
	450	DO j = jts, MIN(jte,jde-1)
	451	DO i = its, MIN(ite,ide-1)
	452	if (LANDMASK(I,J) .gt. 0.5) SM(I,J)=0.
	453	if (LANDMASK(I,J) .le. 0.5) SM(I,J)=1.
	454	if (tsk_gc(I,J) .gt. 0.) then
	455	NMM_TSK(I,J)=tsk_gc(I,J)
	456	else
	457	NMM_TSK(I,J)=t_gc(I,J,1) ! stopgap measure
	458	endif
	459	!
	460	GLAT(I,J)=hlat_gc(I,J)*DEGRAD
	461	GLON(I,J)=hlon_gc(I,J)*DEGRAD
	462	WEASD(I,J)=SNOW(I,J)
	463	XICE(I,J)=XICE_gc(I,J)
	464	ENDDO
	465	ENDDO
	466	! First item is to define the target vertical coordinate
	467
	468	num_metgrid_levels = grid%num_metgrid_levels
	469	eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
	470	ptsgm = model_config_rec%ptsgm
	471	p_top_requested = grid%p_top_requested
	472	pt=p_top_requested
	473
	474	if (internal_time_loop .eq. 1) then
	475
	476	if (eta_levels(1) .ne. 1.0) then
	477	write(message,) '******************************************************************** '
	478	CALL wrf_message(message)
	479	write(message,) '* eta_levels appears not to be specified in the namelist'
	480	CALL wrf_message(message)
	481	write(message,) '* We will call compute_nmm_levels to define layer thicknesses.'
	482	CALL wrf_message(message)
	483	write(message,) '* These levels should be reasonable for running the model, '
	484	CALL wrf_message(message)
	485	write(message,) '* but may not be ideal for the simulation being made. Consider '
	486	CALL wrf_message(message)
	487	write(message,) '* defining your own levels by specifying eta_levels in the model '
	488	CALL wrf_message(message)
	489	write(message,) '* namelist. '
	490	CALL wrf_message(message)
	491	write(message,) '********************************************************************* '
	492	CALL wrf_message(message)
	493
	494	CALL compute_nmm_levels(KDE,p_top_requested,eta_levels)
	495
	496	DO L=1,KDE
	497	write(message,*) 'L, eta_levels(L) returned :: ', L,eta_levels(L)
	498	CALL wrf_message(message)
	499	ENDDO
	500
	501	endif
	502
	503	write(message,*) 'KDE-1: ', KDE-1
	504	CALL wrf_debug(1,message)
	505	allocate(SG1(1:KDE-1))
	506	allocate(SG2(1:KDE-1))
	507	allocate(DSG1(1:KDE-1))
	508	allocate(DSG2(1:KDE-1))
	509	allocate(SGML1(1:KDE))
	510	allocate(SGML2(1:KDE))
	511
	512	CALL define_nmm_vertical_coord (kde-1, ptsgm, pt,pdtop, eta_levels, &
	513	ETA1,DETA1,AETA1, &
	514	ETA2,DETA2,AETA2, DFL, DFRLG )
	515
	516	DO L=KDS,KDE-1
	517	DETA(L)=eta_levels(L)-eta_levels(L+1)
	518	ENDDO
	519	endif
	520
	521	if (.NOT. allocated(PDVP)) allocate(PDVP(IMS:IME,JMS:JME))
	522	if (.NOT. allocated(P3D_OUT)) allocate(P3D_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
	523	if (.NOT. allocated(PSFC_OUTV)) allocate(PSFC_OUTV(IMS:IME,JMS:JME))
	524	if (.NOT. allocated(P3DV_OUT)) allocate(P3DV_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
	525	if (.NOT. allocated(P3DV_IN)) allocate(P3DV_IN(IMS:IME,JMS:JME,num_metgrid_levels))
	526
	527	write(message,*) 'num_metgrid_levels: ', num_metgrid_levels
	528	CALL wrf_message(message)
	529
	530	DO j = jts, MIN(jte,jde-1)
	531	DO i = its, MIN(ite,ide-1)
	532	FIS(I,J)=ht_gc(I,J)*g
	533	!
	534	! IF ( p_gc(I,J,1) .ne. 200100. .AND. (ht_gc(I,J) .eq. ght_gc(I,J,1)) .AND. ht_gc(I,J) .ne. 0) THEN
	535	IF ( p_gc(I,J,1) .ne. 200100. .AND. (abs(ht_gc(I,J)-ght_gc(I,J,1)) .lt. 0.01) .AND. ht_gc(I,J) .ne. 0) THEN
	536	IF (mod(I,10) .eq. 0 .and. mod(J,10) .eq. 0) THEN
	537	write(message,*) 'ht_gc and toposoil to swap, flag_soilhgt ::: ', &
	538	I,J, ht_gc(I,J),toposoil(I,J),flag_soilhgt
	539	CALL wrf_debug(10,message)
	540	ENDIF
	541	IF ( ( flag_soilhgt.EQ. 1 ) ) THEN
	542	ght_gc(I,J,1)=toposoil(I,J)
	543	ENDIF
	544	ENDIF
	545
	546	ENDDO
	547	ENDDO
	548
	549	CALL compute_nmm_surfacep (ht_gc, ght_gc, p_gc , t_gc &
	550	&, psfc_out, num_metgrid_levels &
	551	&, IDS,IDE,JDS,JDE,KDS,KDE &
	552	&, IMS,IME,JMS,JME,KMS,KME &
	553	&, ITS,ITE,JTS,JTE,KTS,KTE ) ! H points
	554
	555	CALL compute_3d_pressure (psfc_out,AETA1,AETA2 &
	556	&, pdtop,pt,pd,p3d_out &
	557	&, IDS,IDE,JDS,JDE,KDS,KDE &
	558	&, IMS,IME,JMS,JME,KMS,KME &
	559	&, ITS,ITE,JTS,JTE,KTS,KTE )
	560
	561	#ifdef DM_PARALLEL
	562	ips=its ; ipe=ite ; jps=jts ; jpe=jte ; kps=kts ; kpe=kte
	563	# include "HALO_NMM_MG2.inc"
	564	#endif
	565
	566	#ifdef DM_PARALLEL
	567	# include "HALO_NMM_MG3.inc"
	568	#endif
	569
	570	do K=1,num_metgrid_levels
	571	do J=JTS,min(JTE,JDE-1)
	572	do I=ITS,min(ITE,IDE-1)
	573
	574	IF (K .eq. KTS) THEN
	575	IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
	576	PDVP(I,J)=0.5*(PD(I,J)+PD(I+1,J))
	577	PSFC_OUTV(I,J)=0.5*(PSFC_OUT(I,J)+PSFC_OUT(I+1,J))
	578	ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
	579	PDVP(I,J)=0.5*(PD(I,J)+PD(I+1,J))
	580	PSFC_OUTV(I,J)=0.5*(PSFC_OUT(I,J)+PSFC_OUT(I+1,J))
	581	ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
	582	PDVP(I,J)=0.5*(PD(I,J-1)+PD(I,J+1))
	583	PSFC_OUTV(I,J)=0.5*(PSFC_OUT(I,J-1)+PSFC_OUT(I,J+1))
	584	ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
	585	PDVP(I,J)=0.5*(PD(I,J-1)+PD(I,J+1))
	586	PSFC_OUTV(I,J)=0.5*(PSFC_OUT(I,J-1)+PSFC_OUT(I,J+1))
	587	ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
	588	PDVP(I,J)=PD(I,J)
	589	PSFC_OUTV(I,J)=PSFC_OUT(I,J)
	590	ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
	591	PDVP(I,J)=0.25*(PD(I,J)+PD(I-1,J)+PD(I,J+1)+PD(I,J-1))
	592	PSFC_OUTV(I,J)=0.25*(PSFC_OUT(I,J)+PSFC_OUT(I-1,J)+ &
	593	PSFC_OUT(I,J+1)+PSFC_OUT(I,J-1))
	594	ELSE ! interior odd row
	595	PDVP(I,J)=0.25*(PD(I,J)+PD(I+1,J)+PD(I,J+1)+PD(I,J-1))
	596	PSFC_OUTV(I,J)=0.25*(PSFC_OUT(I,J)+PSFC_OUT(I+1,J)+ &
	597	PSFC_OUT(I,J+1)+PSFC_OUT(I,J-1))
	598	ENDIF
	599	ENDIF
	600
	601	IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
	602	P3DV_IN(I,J,K)=0.5*(p_gc(I,J,K)+p_gc(I+1,J,K))
	603	ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
	604	P3DV_IN(I,J,K)=0.5*(p_gc(I,J,K)+p_gc(I+1,J,K))
	605	ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
	606	P3DV_IN(I,J,K)=0.5*(p_gc(I,J-1,K)+p_gc(I,J+1,K))
	607	ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
	608	P3DV_IN(I,J,K)=0.5*(p_gc(I,J-1,K)+p_gc(I,J+1,K))
	609	ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
	610	P3DV_IN(I,J,K)=p_gc(I,J,K)
	611	ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
	612	P3DV_IN(I,J,K)=0.25*(p_gc(I,J,K)+p_gc(I-1,J,K) + &
	613	p_gc(I,J+1,K)+p_gc(I,J-1,K))
	614	ELSE ! interior odd row
	615	P3DV_IN(I,J,K)=0.25*(p_gc(I,J,K)+p_gc(I+1,J,K) + &
	616	p_gc(I,J+1,K)+p_gc(I,J-1,K))
	617	ENDIF
	618
	619	enddo
	620	enddo
	621	enddo
	622
	623	CALL compute_3d_pressure (psfc_outv,AETA1,AETA2 &
	624	&, pdtop,pt,pdvp,p3dv_out &
	625	&, IDS,IDE,JDS,JDE,KDS,KDE &
	626	&, IMS,IME,JMS,JME,KMS,KME &
	627	&, ITS,ITE,JTS,JTE,KTS,KTE )
	628
	629	CALL interp_press2press_lin(p_gc, p3d_out &
	630	&, t_gc, T,num_metgrid_levels &
	631	&, .TRUE.,.TRUE.,.TRUE. & ! extrap, ignore_lowest, t_field
	632	&, IDS,IDE,JDS,JDE,KDS,KDE &
	633	&, IMS,IME,JMS,JME,KMS,KME &
	634	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
	635
	636
	637	CALL interp_press2press_lin(p3dv_in, p3dv_out &
	638	&, u_gc, U,num_metgrid_levels &
	639	&, .FALSE.,.TRUE.,.FALSE. &
	640	&, IDS,IDE,JDS,JDE,KDS,KDE &
	641	&, IMS,IME,JMS,JME,KMS,KME &
	642	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
	643
	644	CALL interp_press2press_lin(p3dv_in, p3dv_out &
	645	&, V_gc, V,num_metgrid_levels &
	646	&, .FALSE.,.TRUE.,.FALSE. &
	647	&, IDS,IDE,JDS,JDE,KDS,KDE &
	648	&, IMS,IME,JMS,JME,KMS,KME &
	649	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
	650
	651	IF (hyb_coor) THEN
	652	CALL wind_adjust(p3dv_in,p3dv_out,U_gc,V_gc,U,V &
	653	&, num_metgrid_levels,5000. &
	654	&, IDS,IDE,JDS,JDE,KDS,KDE &
	655	&, IMS,IME,JMS,JME,KMS,KME &
	656	&, ITS,ITE,JTS,JTE,KTS,KTE )
	657	ENDIF
	658
	659
	660	ALLOCATE(qtmp(IMS:IME,JMS:JME,num_metgrid_levels))
	661	ALLOCATE(qtmp2(IMS:IME,JMS:JME,num_metgrid_levels))
	662
	663	CALL rh_to_mxrat (rh_gc, t_gc, p_gc, qtmp , .TRUE. , &
	664	ids , ide , jds , jde , 1 , num_metgrid_levels , &
	665	ims , ime , jms , jme , 1 , num_metgrid_levels , &
	666	its , ite , jts , jte , 1 , num_metgrid_levels )
	667
	668	do K=1,num_metgrid_levels
	669	do J=JTS,min(JTE,JDE-1)
	670	do I=ITS,min(ITE,IDE-1)
	671	QTMP2(I,J,K)=QTMP(I,J,K)/(1.0+QTMP(I,J,K))
	672	end do
	673	end do
	674	end do
	675
	676	CALL interp_press2press_log(p_gc, p3d_out &
	677	&, QTMP2, Q,num_metgrid_levels &
	678	&, .FALSE.,.TRUE. &
	679	&, IDS,IDE,JDS,JDE,KDS,KDE &
	680	&, IMS,IME,JMS,JME,KMS,KME &
	681	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
	682
	683	IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP)
	684	IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP2)
	685
	686	! Get the monthly values interpolated to the current date
	687	! for the traditional monthly
	688	! fields of green-ness fraction and background albedo.
	689
	690	if (internal_time_loop .eq. 1) then
	691
	692	CALL monthly_interp_to_date ( greenfrac_gc , current_date , vegfra , &
	693	ids , ide , jds , jde , kds , kde , &
	694	ims , ime , jms , jme , kms , kme , &
	695	its , ite , jts , jte , kts , kte )
	696
	697	CALL monthly_interp_to_date ( albedo12m_gc , current_date , albbck , &
	698	ids , ide , jds , jde , kds , kde , &
	699	ims , ime , jms , jme , kms , kme , &
	700	its , ite , jts , jte , kts , kte )
	701
	702	! Get the min/max of each i,j for the monthly green-ness fraction.
	703
	704	CALL monthly_min_max ( greenfrac_gc , shdmin , shdmax , &
	705	ids , ide , jds , jde , kds , kde , &
	706	ims , ime , jms , jme , kms , kme , &
	707	its , ite , jts , jte , kts , kte )
	708
	709	! The model expects the green-ness values in percent, not fraction.
	710
	711	DO j = jts, MIN(jte,jde-1)
	712	DO i = its, MIN(ite,ide-1)
	713	!! vegfra(i,j) = vegfra(i,j) * 100.
	714	shdmax(i,j) = shdmax(i,j) * 100.
	715	shdmin(i,j) = shdmin(i,j) * 100.
	716	VEGFRC(I,J)=VEGFRA(I,J)
	717	END DO
	718	END DO
	719
	720	! The model expects the albedo fields as
	721	! a fraction, not a percent. Set the water values to 8%.
	722
	723	DO j = jts, MIN(jte,jde-1)
	724	DO i = its, MIN(ite,ide-1)
	725	if (albbck(i,j) .lt. 5.) then
	726	write(message,*) 'reset albedo to 8%... I,J,albbck:: ', I,J,albbck(I,J)
	727	CALL wrf_debug(10,message)
	728	albbck(I,J)=8.
	729	endif
	730	albbck(i,j) = albbck(i,j) / 100.
	731	snoalb(i,j) = snoalb(i,j) / 100.
	732	IF ( landmask(i,j) .LT. 0.5 ) THEN
	733	albbck(i,j) = 0.08
	734	snoalb(i,j) = 0.08
	735	END IF
	736	albase(i,j)=albbck(i,j)
	737	mxsnal(i,j)=snoalb(i,j)
	738	END DO
	739	END DO
	740
	741	endif
	742
	743	! new deallocs
	744	DEALLOCATE(p3d_out,p3dv_out,p3dv_in)
	745
	746	END IF ! <----- END OF VERTICAL INTERPOLATION PART ---->
	747
	748	if (internal_time_loop .eq. 1) then
	749
	750	!!! WEASD has "snow water equivalent" in mm
	751
	752	DO j = jts, MIN(jte,jde-1)
	753	DO i = its, MIN(ite,ide-1)
	754
	755	IF(SM(I,J).GT.0.9) THEN
	756
	757	IF (XICE(I,J) .gt. 0) then
	758	SI(I,J)=1.0
	759	ENDIF
	760
	761	! SEA
	762	EPSR(I,J)=.97
	763	EMBCK(I,J)=.97
	764	GFFC(I,J)=0.
	765	ALBEDO(I,J)=.06
	766	ALBASE(I,J)=.06
	767	IF(SI (I,J).GT.0. ) THEN
	768	! SEA-ICE
	769	SM(I,J)=0.
	770	SI(I,J)=0.
	771	SICE(I,J)=1.
	772	GFFC(I,J)=0. ! just leave zero as irrelevant
	773	ALBEDO(I,J)=.60
	774	ALBASE(I,J)=.60
	775	ENDIF
	776	ELSE
	777
	778	SI(I,J)=5.0*WEASD(I,J)/1000.
	779	! LAND
	780	EPSR(I,J)=1.0
	781	EMBCK(I,J)=1.0
	782	GFFC(I,J)=0.0 ! just leave zero as irrelevant
	783	SICE(I,J)=0.
	784	SNO(I,J)=SI(I,J)*.20
	785	ENDIF
	786	ENDDO
	787	ENDDO
	788
	789	! DETERMINE ALBEDO OVER LAND
	790	DO j = jts, MIN(jte,jde-1)
	791	DO i = its, MIN(ite,ide-1)
	792	IF(SM(I,J).LT.0.9.AND.SICE(I,J).LT.0.9) THEN
	793	! SNOWFREE ALBEDO
	794	IF ( (SNO(I,J) .EQ. 0.0) .OR. &
	795	(ALBASE(I,J) .GE. MXSNAL(I,J) ) ) THEN
	796	ALBEDO(I,J) = ALBASE(I,J)
	797	ELSE
	798	! MODIFY ALBEDO IF SNOWCOVER:
	799	! BELOW SNOWDEPTH THRESHOLD...
	800	IF (SNO(I,J) .LT. SNUP) THEN
	801	RSNOW = SNO(I,J)/SNUP
	802	SNOFAC = 1. - ( EXP(-SALPRSNOW) - RSNOWEXP(-SALP))
	803	! ABOVE SNOWDEPTH THRESHOLD...
	804	ELSE
	805	SNOFAC = 1.0
	806	ENDIF
	807	! CALCULATE ALBEDO ACCOUNTING FOR SNOWDEPTH AND VGFRCK
	808	ALBEDO(I,J) = ALBASE(I,J) &
	809	+ (1.0-VEGFRA(I,J))SNOFAC(MXSNAL(I,J)-ALBASE(I,J))
	810	ENDIF
	811	END IF
	812	SI(I,J)=5.0*WEASD(I,J)
	813	SNO(I,J)=WEASD(I,J)
	814
	815	!! convert VEGFRA
	816	VEGFRA(I,J)=VEGFRA(I,J)*100.
	817	!
	818	ENDDO
	819	ENDDO
	820
	821	#ifdef DM_PARALLEL
	822
	823	ALLOCATE(SM_G(IDS:IDE,JDS:JDE),SICE_G(IDS:IDE,JDS:JDE))
	824
	825	CALL WRF_PATCH_TO_GLOBAL_REAL( SICE(IMS,JMS) &
	826	&, SICE_G,grid%DOMDESC &
	827	&, 'z','xy' &
	828	&, IDS,IDE-1,JDS,JDE-1,1,1 &
	829	&, IMS,IME,JMS,JME,1,1 &
	830	&, ITS,ITE,JTS,JTE,1,1 )
	831
	832	CALL WRF_PATCH_TO_GLOBAL_REAL( SM(IMS,JMS) &
	833	&, SM_G,grid%DOMDESC &
	834	&, 'z','xy' &
	835	&, IDS,IDE-1,JDS,JDE-1,1,1 &
	836	&, IMS,IME,JMS,JME,1,1 &
	837	&, ITS,ITE,JTS,JTE,1,1 )
	838
	839
	840	IF (WRF_DM_ON_MONITOR()) THEN
	841
	842	637 format(40(f3.0,1x))
	843
	844	allocate(IHE_G(JDS:JDE-1),IHW_G(JDS:JDE-1))
	845	DO j = JDS, JDE-1
	846	IHE_G(J)=MOD(J+1,2)
	847	IHW_G(J)=IHE_G(J)-1
	848	ENDDO
	849
	850	DO ITER=1,10
	851	DO j = jds+1, (jde-1)-1
	852	DO i = ids+1, (ide-1)-1
	853
	854	! any sea ice around point in question?
	855
	856	IF (SM_G(I,J) .ge. 0.9) THEN
	857	SEAICESUM=SICE_G(I+IHE_G(J),J+1)+SICE_G(I+IHW_G(J),J+1)+ &
	858	SICE_G(I+IHE_G(J),J-1)+SICE_G(I+IHW_G(J),J-1)
	859	IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
	860
	861	IF ((SICE_G(I+IHE_G(J),J+1).lt.0.1 .and. SM_G(I+IHE_G(J),J+1).lt.0.1) .OR. &
	862	(SICE_G(I+IHW_G(J),J+1).lt.0.1 .and. SM_G(I+IHW_G(J),J+1).lt.0.1) .OR. &
	863	(SICE_G(I+IHE_G(J),J-1).lt.0.1 .and. SM_G(I+IHE_G(J),J-1).lt.0.1) .OR. &
	864	(SICE_G(I+IHW_G(J),J-1).lt.0.1 .and. SM_G(I+IHW_G(J),J-1).lt.0.1)) THEN
	865
	866	! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
	867
	868	write(message,*) 'making seaice (1): ', I,J
	869	CALL wrf_debug(100,message)
	870	SICE_G(I,J)=1.0
	871	SM_G(I,J)=0.
	872
	873	ENDIF
	874
	875	ELSEIF (SEAICESUM .ge. 3) THEN
	876
	877	! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
	878
	879	write(message,*) 'making seaice (2): ', I,J
	880	CALL wrf_debug(100,message)
	881	SICE_G(I,J)=1.0
	882	SM_G(I,J)=0.
	883	ENDIF
	884
	885	ENDIF
	886
	887	ENDDO
	888	ENDDO
	889	ENDDO
	890
	891	ENDIF
	892
	893	CALL WRF_GLOBAL_TO_PATCH_REAL( SICE_G, SICE &
	894	&, grid%DOMDESC &
	895	&, 'z','xy' &
	896	&, IDS,IDE-1,JDS,JDE-1,1,1 &
	897	&, IMS,IME,JMS,JME,1,1 &
	898	&, ITS,ITE,JTS,JTE,1,1 )
	899
	900	CALL WRF_GLOBAL_TO_PATCH_REAL( SM_G,SM &
	901	&, grid%DOMDESC &
	902	&, 'z','xy' &
	903	&, IDS,IDE-1,JDS,JDE-1,1,1 &
	904	&, IMS,IME,JMS,JME,1,1 &
	905	&, ITS,ITE,JTS,JTE,1,1 )
	906
	907	IF (WRF_DM_ON_MONITOR()) THEN
	908
	909	DEALLOCATE(SM_G,SICE_G)
	910	DEALLOCATE(IHE_G,IHW_G)
	911
	912	ENDIF
	913
	914	! write(message,*) 'revised sea ice on patch'
	915	! CALL wrf_debug(100,message)
	916	! DO J=JTE,JTS,-(((JTE-JTS)/25)+1)
	917	! write(message,637) (SICE(I,J),I=ITS,ITE,ITE/20)
	918	! CALL wrf_debug(100,message)
	919	! END DO
	920
	921	#else
	922	! serial sea ice reprocessing
	923
	924	DO j = jts, MIN(jte,jde-1)
	925	IHE(J)=MOD(J+1,2)
	926	IHW(J)=IHE(J)-1
	927	ENDDO
	928
	929	DO ITER=1,10
	930	DO j = jts+1, MIN(jte,jde-1)-1
	931	DO i = its+1, MIN(ite,ide-1)-1
	932
	933	! any sea ice around point in question?
	934
	935	IF (SM(I,J) .gt. 0.9) THEN
	936	SEAICESUM=SICE(I+IHE(J),J+1)+SICE(I+IHW(J),J+1)+ &
	937	SICE(I+IHE(J),J-1)+SICE(I+IHW(J),J-1)
	938	IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
	939	IF ((SICE(I+IHE(J),J+1).lt.0.1 .and. SM(I+IHE(J),J+1).lt.0.1) .OR. &
	940	(SICE(I+IHW(J),J+1).lt.0.1 .and. SM(I+IHW(J),J+1).lt.0.1) .OR. &
	941	(SICE(I+IHE(J),J-1).lt.0.1 .and. SM(I+IHE(J),J-1).lt.0.1) .OR. &
	942	(SICE(I+IHW(J),J-1).lt.0.1 .and. SM(I+IHW(J),J-1).lt.0.1)) THEN
	943
	944	! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
	945	SICE(I,J)=1.0
	946	SM(I,J)=0.
	947	ENDIF
	948	ELSEIF (SEAICESUM .ge. 3) THEN
	949	! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
	950	SICE(I,J)=1.0
	951	SM(I,J)=0.
	952	ENDIF
	953	ENDIF
	954
	955	ENDDO
	956	ENDDO
	957	ENDDO
	958
	959	#endif
	960
	961	! this block meant to guarantee land/sea agreement between SM and landmask
	962
	963	DO j = jts, MIN(jte,jde-1)
	964	DO i = its, MIN(ite,ide-1)
	965
	966	IF (SM(I,J) .gt. 0.5) THEN
	967	landmask(I,J)=0.0
	968	ELSEIF (SM(I,J) .lt. 0.5 .and. SICE(I,J) .gt. 0.9) then
	969	landmask(I,J)=0.0
	970	ELSEIF (SM(I,J) .lt. 0.5 .and. SICE(I,J) .lt. 0.1) then
	971	landmask(I,J)=1.0
	972	ELSE
	973	write(message,*) 'missed point in landmask definition ' , I,J
	974	CALL wrf_message(message)
	975	landmask(I,J)=0.0
	976	ENDIF
	977	!
	978	IF (SICE(I,J) .gt. 0.5 .and. NMM_TSK(I,J) .lt. 0.1 .and. SST(I,J) .gt. 0.) THEN
	979	write(message,*) 'set NMM_TSK to: ', SST(I,J)
	980	CALL wrf_message(message)
	981	NMM_TSK(I,J)=SST(I,J)
	982	SST(I,J)=0.
	983	endif
	984
	985	ENDDO
	986	ENDDO
	987
	988	! For sf_surface_physics = 1, we want to use close to a 10 cm value
	989	! for the bottom level of the soil temps.
	990
	991	IF ( ( model_config_rec%sf_surface_physics(grid%id) .EQ. 1 ) .AND. &
	992	( flag_st000010 .EQ. 1 ) ) THEN
	993	DO j = jts , MIN(jde-1,jte)
	994	DO i = its , MIN(ide-1,ite)
	995	soiltb(i,j) = st000010(i,j)
	996	END DO
	997	END DO
	998	END IF
	999
	1000	! Adjust the various soil temperature values depending on the difference in
	1001	! in elevation between the current model's elevation and the incoming data's
	1002	! orography.
	1003
	1004	IF ( ( flag_toposoil .EQ. 1 ) ) THEN
	1005
	1006	ALLOCATE(HT(ims:ime,jms:jme))
	1007
	1008	DO J=jms,jme
	1009	DO I=ims,ime
	1010	HT(I,J)=FIS(I,J)/9.81
	1011	END DO
	1012	END DO
	1013
	1014	! if (maxval(toposoil) .gt. 100.) then
	1015	!
	1016	! Being avoided. Something to revisit eventually.
	1017	!
	1018	!1219 might be simply a matter of including TOPOSOIL
	1019	!
	1020	! CODE NOT TESTED AT NCEP USING THIS FUNCTIONALITY,
	1021	! SO TO BE SAFE WILL AVOID FOR RETRO RUNS.
	1022	!
	1023	! CALL adjust_soil_temp_new ( soiltb , 2 , &
	1024	! nmm_tsk , ht , toposoil , landmask, flag_toposoil , &
	1025	! st000010 , st010040 , st040100 , st100200 , st010200 , &
	1026	! flag_st000010 , flag_st010040 , flag_st040100 , &
	1027	! flag_st100200 , flag_st010200 , &
	1028	! soilt000 , soilt005 , soilt020 , soilt040 , &
	1029	! soilt160 , soilt300 , &
	1030	! flag_soilt000 , flag_soilt005 , flag_soilt020 , &
	1031	! flag_soilt040 , flag_soilt160 , flag_soilt300 , &
	1032	! ids , ide , jds , jde , kds , kde , &
	1033	! ims , ime , jms , jme , kms , kme , &
	1034	! its , ite , jts , jte , kts , kte )
	1035	! endif
	1036
	1037	END IF
	1038
	1039	! Process the LSM data.
	1040
	1041	! surface_input_source=1 => use data from static file
	1042	! (fractional category as input)
	1043	! surface_input_source=2 => use data from grib file
	1044	! (dominant category as input)
	1045
	1046	IF ( config_flags%surface_input_source .EQ. 1 ) THEN
	1047	vegcat (its,jts) = 0
	1048	soilcat(its,jts) = 0
	1049	END IF
	1050
	1051	! Generate the vegetation and soil category information
	1052	! from the fractional input
	1053	! data, or use the existing dominant category fields if they exist.
	1054
	1055	IF ((soilcat(its,jts) .LT. 0.5) .AND. (vegcat(its,jts) .LT. 0.5)) THEN
	1056
	1057	num_veg_cat = SIZE ( landusef_gc , DIM=3 )
	1058	num_soil_top_cat = SIZE ( soilctop_gc , DIM=3 )
	1059	num_soil_bot_cat = SIZE ( soilcbot_gc , DIM=3 )
	1060
	1061	do J=JMS,JME
	1062	do K=1,num_veg_cat
	1063	do I=IMS,IME
	1064	landusef(I,K,J)=landusef_gc(I,J,K)
	1065	enddo
	1066	enddo
	1067	enddo
	1068
	1069	do J=JMS,JME
	1070	do K=1,num_soil_top_cat
	1071	do I=IMS,IME
	1072	soilctop(I,K,J)=soilctop_gc(I,J,K)
	1073	enddo
	1074	enddo
	1075	enddo
	1076
	1077	do J=JMS,JME
	1078	do K=1,num_soil_bot_cat
	1079	do I=IMS,IME
	1080	soilcbot(I,K,J)=soilcbot_gc(I,J,K)
	1081	enddo
	1082	enddo
	1083	enddo
	1084
	1085	! sm (1=water, 0=land)
	1086	! landmask(0=water, 1=land)
	1087
	1088
	1089	write(message,*) 'landmask into process_percent_cat_new'
	1090
	1091	CALL wrf_debug(1,message)
	1092	do J=JTE,JTS,-(((JTE-JTS)/20)+1)
	1093	write(message,641) (landmask(I,J),I=ITS,min(ITE,IDE-1),((ITE-ITS)/15)+1)
	1094	CALL wrf_debug(1,message)
	1095	enddo
	1096	641 format(25(f3.0,1x))
	1097
	1098	CALL process_percent_cat_new ( landmask , &
	1099	landusef , soilctop , soilcbot , &
	1100	isltyp , ivgtyp , &
	1101	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
	1102	ids , ide , jds , jde , kds , kde , &
	1103	ims , ime , jms , jme , kms , kme , &
	1104	its , ite , jts , jte , kts , kte , &
	1105	model_config_rec%iswater(grid%id) )
	1106
	1107	DO j = jts , MIN(jde-1,jte)
	1108	DO i = its , MIN(ide-1,ite)
	1109	vegcat(i,j) = ivgtyp(i,j)
	1110	soilcat(i,j) = isltyp(i,j)
	1111	END DO
	1112	END DO
	1113
	1114	ELSE
	1115
	1116	! Do we have dominant soil and veg data from the input already?
	1117
	1118	IF ( soilcat(its,jts) .GT. 0.5 ) THEN
	1119	DO j = jts, MIN(jde-1,jte)
	1120	DO i = its, MIN(ide-1,ite)
	1121	isltyp(i,j) = NINT( soilcat(i,j) )
	1122	END DO
	1123	END DO
	1124	END IF
	1125	IF ( vegcat(its,jts) .GT. 0.5 ) THEN
	1126	DO j = jts, MIN(jde-1,jte)
	1127	DO i = its, MIN(ide-1,ite)
	1128	ivgtyp(i,j) = NINT( vegcat(i,j) )
	1129	END DO
	1130	END DO
	1131	END IF
	1132
	1133	ENDIF
	1134
	1135	DO j = jts, MIN(jde-1,jte)
	1136	DO i = its, MIN(ide-1,ite)
	1137
	1138	IF (SICE(I,J) .lt. 0.1) THEN
	1139	IF (landmask(I,J) .gt. 0.5 .and. sm(I,J) .gt. 0.5) THEN
	1140	write(message,*) 'land mask and SM both > 0.5: ', &
	1141	I,J,landmask(I,J),sm(I,J)
	1142	CALL wrf_message(message)
	1143	SM(I,J)=0.
	1144	ELSEIF (landmask(I,J) .lt. 0.5 .and. sm(I,J) .lt. 0.5) THEN
	1145	write(message,*) 'land mask and SM both < 0.5: ', &
	1146	I,J, landmask(I,J),sm(I,J)
	1147	CALL wrf_message(message)
	1148	SM(I,J)=1.
	1149	ENDIF
	1150	ELSE
	1151	IF (landmask(I,J) .gt. 0.5 .and. SM(I,J)+SICE(I,J) .gt. 0.9) then
	1152	write(message,*) 'landmask says LAND, SM/SICE say SEAICE: ', I,J
	1153	ENDIF
	1154	ENDIF
	1155
	1156	ENDDO
	1157	ENDDO
	1158
	1159	DO j = jts, MIN(jde-1,jte)
	1160	DO i = its, MIN(ide-1,ite)
	1161
	1162	if (SICE(I,J) .gt. 0.9) then
	1163	ISLTYP(I,J)=16
	1164	IVGTYP(I,J)=24
	1165	endif
	1166
	1167	ENDDO
	1168	ENDDO
	1169
	1170	DO j = jts, MIN(jde-1,jte)
	1171	DO i = its, MIN(ide-1,ite)
	1172
	1173	if (SM(I,J) .lt. 0.5) then
	1174	SST(I,J)=0.
	1175	endif
	1176
	1177	if (SM(I,J) .gt. 0.5) then
	1178	if (SST(I,J) .lt. 0.1) then
	1179	SST(I,J)=NMM_TSK(I,J)
	1180	endif
	1181	NMM_TSK(I,J)=0.
	1182	endif
	1183
	1184	IF ( (NMM_TSK(I,J)+SST(I,J)) .lt. 200. .or. &
	1185	(NMM_TSK(I,J)+SST(I,J)) .gt. 350. ) THEN
	1186	write(message,*) 'TSK, SST trouble at : ', I,J
	1187	CALL wrf_message(message)
	1188	write(message,*) 'SM, NMM_TSK,SST ', SM(I,J),NMM_TSK(I,J),SST(I,J)
	1189	CALL wrf_message(message)
	1190	ENDIF
	1191
	1192	ENDDO
	1193	ENDDO
	1194
	1195	write(message,*) 'SM'
	1196	CALL wrf_message(message)
	1197
	1198	DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
	1199	write(message,635) (sm(i,J),I=its,ite,((ite-its)/10)+1)
	1200	CALL wrf_message(message)
	1201	END DO
	1202
	1203	write(message,*) 'SST/NMM_TSK'
	1204	CALL wrf_debug(10,message)
	1205	DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
	1206	write(message,635) (SST(I,J)+NMM_TSK(I,J),I=ITS,min(ide-1,ite),((ite-its)/10)+1)
	1207	CALL wrf_debug(10,message)
	1208	END DO
	1209
	1210	635 format(20(f5.1,1x))
	1211
	1212	DO j = jts, MIN(jde-1,jte)
	1213	DO i = its, MIN(ide-1,ite)
	1214	IF ( ( landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
	1215	soiltb(i,j) = sst(i,j)
	1216	ELSE IF ( landmask(i,j) .GT. 0.5 ) THEN
	1217	soiltb(i,j) = nmm_tsk(i,j)
	1218	END IF
	1219	END DO
	1220	END DO
	1221
	1222	! END IF
	1223
	1224	! Land use categories, dominant soil and vegetation types (if available).
	1225
	1226	! allocate(lu_index(ims:ime,jms:jme))
	1227
	1228	DO j = jts, MIN(jde-1,jte)
	1229	DO i = its, MIN(ide-1,ite)
	1230	lu_index(i,j) = ivgtyp(i,j)
	1231	END DO
	1232	END DO
	1233
	1234	if (flag_sst .eq. 1) log_flag_sst=.true.
	1235	if (flag_sst .eq. 0) log_flag_sst=.false.
	1236
	1237	write(message,*) 'st_input dimensions: ', size(st_input,dim=1), &
	1238	size(st_input,dim=2),size(st_input,dim=3)
	1239	CALL wrf_debug(100,message)
	1240
	1241	! write(message,*) 'maxval st_input(1): ', maxval(st_input(:,1,:))
	1242	! CALL wrf_message(message)
	1243	! write(message,*) 'maxval st_input(2): ', maxval(st_input(:,2,:))
	1244	! CALL wrf_message(message)
	1245	! write(message,*) 'maxval st_input(3): ', maxval(st_input(:,3,:))
	1246	! CALL wrf_message(message)
	1247	! write(message,*) 'maxval st_input(4): ', maxval(st_input(:,4,:))
	1248	! CALL wrf_message(message)
	1249
	1250	! =============================================================
	1251
	1252	IF (.NOT. ALLOCATED(TG_ALT))ALLOCATE(TG_ALT(grid%sm31:grid%em31,jms:jme))
	1253
	1254	TPH0=TPH0D*DTR
	1255	WBD=-(((ide-1)-1)*DLMD)
	1256	WB= WBD*DTR
	1257	SBD=-(((jde-1)/2)*DPHD)
	1258	SB= SBD*DTR
	1259	DLM=DLMD*DTR
	1260	DPH=DPHD*DTR
	1261	TDLM=DLM+DLM
	1262	TDPH=DPH+DPH
	1263	WBI=WB+TDLM
	1264	SBI=SB+TDPH
	1265	EBI=WB+(ide-2)*TDLM
	1266	ANBI=SB+(jde-2)*DPH
	1267	STPH0=SIN(TPH0)
	1268	CTPH0=COS(TPH0)
	1269	TSPH=3600./GRID%DT
	1270	DO J=JTS,min(JTE,JDE-1)
	1271	TLM=WB-TDLM+MOD(J,2)*DLM !For velocity points on the E grid
	1272	TPH=SB+float(J-1)*DPH
	1273	STPH=SIN(TPH)
	1274	CTPH=COS(TPH)
	1275	DO I=ITS,MIN(ITE,IDE-1)
	1276
	1277	if (I .eq. ITS) THEN
	1278	TLM=TLM+TDLM*ITS
	1279	else
	1280	TLM=TLM+TDLM
	1281	endif
	1282
	1283	TERM1=(STPH0CTPHCOS(TLM)+CTPH0*STPH)
	1284	FP=TWOM*(TERM1)
	1285	F(I,J)=0.5GRID%DTFP
	1286	ENDDO
	1287	ENDDO
	1288	DO J=JTS,min(JTE,JDE-1)
	1289	TLM=WB-TDLM+MOD(J+1,2)*DLM !For mass points on the E grid
	1290	TPH=SB+float(J-1)*DPH
	1291	STPH=SIN(TPH)
	1292	CTPH=COS(TPH)
	1293	DO I=ITS,MIN(ITE,IDE-1)
	1294
	1295	if (I .eq. ITS) THEN
	1296	TLM=TLM+TDLM*ITS
	1297	else
	1298	TLM=TLM+TDLM
	1299	endif
	1300
	1301	TERM1=(STPH0CTPHCOS(TLM)+CTPH0*STPH)
	1302	TERM1=MIN(TERM1,1.0D0)
	1303	TERM1=MAX(TERM1,-1.0D0)
	1304	APH=ASIN(TERM1)
	1305	TG_ALT(I,J)=TG0+TGA*COS(APH)-FIS(I,J)/3333.
	1306	ENDDO
	1307	ENDDO
	1308
	1309	DO j = jts, MIN(jde-1,jte)
	1310	DO i = its, MIN(ide-1,ite)
	1311	! IF ( ( landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
	1312	! SICE(I,J) .eq. 0. ) THEN
	1313	! TG(i,j) = sst(i,j)
	1314	! ELSEIF (SICE(I,J) .eq. 1) THEN
	1315	! TG(i,j) = 271.16
	1316	! END IF
	1317
	1318	if (TG(I,J) .lt. 200.) then ! only use default TG_ALT definition if
	1319	! not getting TGROUND from SI
	1320	TG(I,J)=TG_ALT(I,J)
	1321	endif
	1322
	1323	if (TG(I,J) .lt. 200. .or. TG(I,J) .gt. 320.) then
	1324	write(message,*) 'problematic TG point at : ', I,J
	1325	CALL wrf_message( message )
	1326	endif
	1327
	1328	adum2d(i,j)=nmm_tsk(I,J)+sst(I,J)
	1329
	1330	END DO
	1331	END DO
	1332
	1333	DEALLOCATE(TG_ALT)
	1334
	1335	write(message,*) 'call process_soil_real with num_st_levels_input: ', num_st_levels_input
	1336	CALL wrf_message( message )
	1337
	1338	! =============================================================
	1339
	1340	CALL process_soil_real ( adum2d, TG , &
	1341	landmask, sst, &
	1342	st_input, sm_input, sw_input, &
	1343	st_levels_input , sm_levels_input , &
	1344	sw_levels_input , &
	1345	sldpth , dzsoil , stc , smc , sh2o, &
	1346	flag_sst , flag_soilt000, flag_soilm000, &
	1347	ids , ide , jds , jde , kds , kde , &
	1348	ims , ime , jms , jme , kms , kme , &
	1349	its , ite , jts , jte , kts , kte , &
	1350	model_config_rec%sf_surface_physics(grid%id) , &
	1351	model_config_rec%num_soil_layers , &
	1352	model_config_rec%real_data_init_type , &
	1353	num_st_levels_input , num_sm_levels_input , &
	1354	num_sw_levels_input , &
	1355	num_st_levels_alloc , num_sm_levels_alloc , &
	1356	num_sw_levels_alloc )
	1357
	1358	! =============================================================
	1359
	1360	! Minimum soil values, residual, from RUC LSM scheme.
	1361	! For input from Noah and using
	1362	! RUC LSM scheme, this must be subtracted from the input
	1363	! total soil moisture. For input RUC data and using the Noah LSM scheme,
	1364	! this value must be added to the soil moisture_input.
	1365
	1366	lqmi(1:num_soil_top_cat) = &
	1367	(/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
	1368	0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
	1369	0.004, 0.065 /) !dusan , 0.020, 0.004, 0.008 /)
	1370
	1371	! At the initial time we care about values of soil moisture and temperature,
	1372	! other times are ignored by the model, so we ignore them, too.
	1373
	1374	account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
	1375
	1376	CASE ( LSMSCHEME , NMMLSMSCHEME)
	1377	iicount = 0
	1378	IF ( FLAG_SM000010 .EQ. 1 ) THEN
	1379	DO j = jts, MIN(jde-1,jte)
	1380	DO i = its, MIN(ide-1,ite)
	1381	IF ((landmask(i,j).gt.0.5) .and. (stc(i,1,j) .gt. 200) .and. &
	1382	(stc(i,1,j) .lt. 400) .and. (smc(i,1,j) .lt. 0.005)) then
	1383	write(message,*) 'Noah > Noah: bad soil moisture at i,j = ',i,j,smc(i,:,j)
	1384	CALL wrf_message(message)
	1385	iicount = iicount + 1
	1386	smc(i,:,j) = 0.005
	1387	END IF
	1388	END DO
	1389	END DO
	1390	IF ( iicount .GT. 0 ) THEN
	1391	write(message,*) 'Noah -> Noah: total number of small soil moisture locations= ',&
	1392	iicount
	1393	CALL wrf_message(message)
	1394	END IF
	1395	ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
	1396	DO j = jts, MIN(jde-1,jte)
	1397	DO i = its, MIN(ide-1,ite)
	1398	smc(i,:,j) = smc(i,:,j) + lqmi(isltyp(i,j))
	1399	END DO
	1400	END DO
	1401	DO j = jts, MIN(jde-1,jte)
	1402	DO i = its, MIN(ide-1,ite)
	1403	IF ((landmask(i,j).gt.0.5) .and. (stc(i,1,j) .gt. 200) .and. &
	1404	(stc(i,1,j) .lt. 400) .and. (smc(i,1,j) .lt. 0.004)) then
	1405	write(message,*) 'RUC -> Noah: bad soil moisture at i,j = ' &
	1406	,i,j,smc(i,:,j)
	1407	CALL wrf_message(message)
	1408	iicount = iicount + 1
	1409	smc(i,:,j) = 0.004
	1410	END IF
	1411	END DO
	1412	END DO
	1413	IF ( iicount .GT. 0 ) THEN
	1414	write(message,*) 'RUC -> Noah: total number of small soil moisture locations = ',&
	1415	iicount
	1416	CALL wrf_message(message)
	1417	END IF
	1418	END IF
	1419	CASE ( RUCLSMSCHEME )
	1420	iicount = 0
	1421	IF ( FLAG_SM000010 .EQ. 1 ) THEN
	1422	DO j = jts, MIN(jde-1,jte)
	1423	DO i = its, MIN(ide-1,ite)
	1424	smc(i,:,j) = MAX ( smc(i,:,j) - lqmi(isltyp(i,j)) , 0. )
	1425	END DO
	1426	END DO
	1427	ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
	1428	! no op
	1429	END IF
	1430
	1431	END SELECT account_for_zero_soil_moisture
	1432
	1433	!!! zero out NMM_TSK at water points again
	1434
	1435	DO j = jts, MIN(jde-1,jte)
	1436	DO i = its, MIN(ide-1,ite)
	1437	if (SM(I,J) .gt. 0.5) then
	1438	NMM_TSK(I,J)=0.
	1439	endif
	1440	END DO
	1441	END DO
	1442
	1443	!! check on STC
	1444
	1445	DO j = jts, MIN(jde-1,jte)
	1446	DO i = its, MIN(ide-1,ite)
	1447
	1448	IF (SICE(I,J) .gt. 0.9) then
	1449	DO L = 1, grid%num_soil_layers
	1450	STC(I,L,J)=271.16 ! TG value used by Eta/NMM
	1451	END DO
	1452	END IF
	1453
	1454	IF (SM(I,J) .gt. 0.9) then
	1455	DO L = 1, grid%num_soil_layers
	1456	STC(I,L,J)=273.16 ! TG value used by Eta/NMM
	1457	END DO
	1458	END IF
	1459
	1460	END DO
	1461	END DO
	1462
	1463	DO j = jts, MIN(jde-1,jte)
	1464	DO i = its, MIN(ide-1,ite)
	1465
	1466	if (SM(I,J) .lt. 0.1 .and. STC(I,1,J) .lt. 0.1) THEN
	1467	write(message,*) 'troublesome SM,STC,SMC value: ', I,J,SM(I,J), stc(I,1,J),smc(I,1,J)
	1468	CALL wrf_message(message)
	1469	do JJ=J-1,J+1
	1470	do L=1, grid%num_soil_layers
	1471	do II=I-1,I+1
	1472
	1473	if (II .ge. its .and. II .le. MIN(ide-1,ite) .and. &
	1474	JJ .ge. jts .and. JJ .le. MIN(jde-1,jte)) then
	1475
	1476	STC(I,L,J)=amax1(STC(I,L,J),STC(II,L,JJ))
	1477	cur_smc=SMC(I,L,J)
	1478
	1479	if ( SMC(II,L,JJ) .gt. 0.005 .and. SMC(II,L,JJ) .lt. 1.0) then
	1480	aposs_smc=SMC(II,L,JJ)
	1481
	1482	if ( cur_smc .eq. 0 ) then
	1483	cur_smc=aposs_smc
	1484	SMC(I,L,J)=cur_smc
	1485	else
	1486	cur_smc=amin1(cur_smc,aposs_smc)
	1487	cur_smc=amin1(cur_smc,aposs_smc)
	1488	SMC(I,L,J)=cur_smc
	1489	endif
	1490	endif
	1491
	1492	endif ! bounds check
	1493
	1494	enddo
	1495	enddo
	1496	enddo
	1497	write(message,*) 'STC, SMC(1) now: ', stc(I,1,J),smc(I,1,J)
	1498	CALL wrf_message(message)
	1499	endif
	1500
	1501	if (STC(I,1,J) .lt. 0.1) then
	1502	write(message,*) 'QUITTING DUE TO STILL troublesome STC value: ', I,J, stc(I,1,J),smc(I,1,J)
	1503	call wrf_error_fatal(message)
	1504	endif
	1505
	1506	ENDDO
	1507	ENDDO
	1508
	1509	!hardwire soil stuff for time being
	1510
	1511	! RTDPTH=0.
	1512	! RTDPTH(1)=0.1
	1513	! RTDPTH(2)=0.3
	1514	! RTDPTH(3)=0.6
	1515
	1516	! SLDPTH=0.
	1517	! SLDPTH(1)=0.1
	1518	! SLDPTH(2)=0.3
	1519	! SLDPTH(3)=0.6
	1520	! SLDPTH(4)=1.0
	1521
	1522	!!! main body of nmm_specific starts here
	1523	!
	1524	do J=jts,min(jte,jde-1)
	1525	do I=its,min(ite,ide-1)
	1526	RES(I,J)=1.
	1527	enddo
	1528	enddo
	1529
	1530	!! HBM2
	1531
	1532	HBM2=0.
	1533
	1534	do J=jts,min(jte,jde-1)
	1535	do I=its,min(ite,ide-1)
	1536
	1537	IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
	1538	(I .ge. 2 .and. I .le. (ide-1)-2+mod(J,2)) ) THEN
	1539	HBM2(I,J)=1.
	1540	ENDIF
	1541	enddo
	1542	enddo
	1543
	1544	!! HBM3
	1545	HBM3=0.
	1546
	1547	!! LOOP OVER LOCAL DIMENSIONS
	1548
	1549	do J=jts,min(jte,jde-1)
	1550	IHWG(J)=mod(J+1,2)-1
	1551	IF (J .ge. 4 .and. J .le. (jde-1)-3) THEN
	1552	IHL=(ids+1)-IHWG(J)
	1553	IHH=(ide-1)-2
	1554	do I=its,min(ite,ide-1)
	1555	IF (I .ge. IHL .and. I .le. IHH) HBM3(I,J)=1.
	1556	enddo
	1557	ENDIF
	1558	enddo
	1559
	1560	!! VBM2
	1561
	1562	VBM2=0.
	1563
	1564	do J=jts,min(jte,jde-1)
	1565	do I=its,min(ite,ide-1)
	1566
	1567	IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
	1568	(I .ge. 2 .and. I .le. (ide-1)-1-mod(J,2)) ) THEN
	1569
	1570	VBM2(I,J)=1.
	1571
	1572	ENDIF
	1573
	1574	enddo
	1575	enddo
	1576
	1577	!! VBM3
	1578
	1579	VBM3=0.
	1580
	1581	do J=jts,min(jte,jde-1)
	1582	do I=its,min(ite,ide-1)
	1583
	1584	IF ( (J .ge. 4 .and. J .le. (jde-1)-3) .AND. &
	1585	(I .ge. 3-mod(J,2) .and. I .le. (ide-1)-2) ) THEN
	1586	VBM3(I,J)=1.
	1587	ENDIF
	1588
	1589	enddo
	1590	enddo
	1591
	1592	DTAD=1.0
	1593	! IDTCF=DTCF, IDTCF=4
	1594	DTCF=4.0 ! used?
	1595
	1596	DY_NMM=ERAD*DPH
	1597	CPGFV=-GRID%DT/(48.*DY_NMM)
	1598	EN= GRID%DT/( 4.DY_NMM)DTAD
	1599	ENT=GRID%DT/(16.DY_NMM)DTAD
	1600
	1601	DO J=jts,nnyp
	1602	KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2
	1603	KVL2(J)=(IDE-1)*(J-1)-J/2+2
	1604	KHH2(J)=(IDE-1)*J-J/2-1
	1605	KVH2(J)=(IDE-1)*J-(J+1)/2-1
	1606	ENDDO
	1607
	1608	TPH=SB-DPH
	1609
	1610	DO J=jts,min(jte,jde-1)
	1611	TPH=SB+float(J-1)*DPH
	1612	DXP=ERADDLMCOS(TPH)
	1613	DXJ(J)=DXP
	1614	WPDARJ(J)=-W_NMM * &
	1615	((ERADDLMAMIN1(COS(ANBI),COS(SBI)))2+DY_NMM2)/ &
	1616	(GRID%DT32.DXP*DY_NMM)
	1617
	1618	CPGFUJ(J)=-GRID%DT/(48.*DXP)
	1619	CURVJ(J)=.5GRID%DTTAN(TPH)/ERAD
	1620	FCPJ(J)=GRID%DT/(CP192.DXP*DY_NMM)
	1621	FDIVJ(J)=1./(12.DXPDY_NMM)
	1622	! EMJ(J)= GRID%DT/( 4.DXP)DTAD
	1623	! EMTJ(J)=GRID%DT/(16.DXP)DTAD
	1624	FADJ(J)=-GRID%DT/(48.DXPDY_NMM)*DTAD
	1625	ACDT=GRID%DTSQRT((ERADDLMAMIN1(COS(ANBI),COS(SBI)))2+DY_NMM*2)
	1626	CDDAMP=CODAMP*ACDT
	1627	HDACJ(J)=COACACDT/(4.DXP*DY_NMM)
	1628	DDMPUJ(J)=CDDAMP/DXP
	1629	DDMPVJ(J)=CDDAMP/DY_NMM
	1630	ENDDO
	1631
	1632	DO J=JTS,min(JTE,JDE-1)
	1633	TLM=WB-TDLM+MOD(J,2)*DLM
	1634	TPH=SB+float(J-1)*DPH
	1635	STPH=SIN(TPH)
	1636	CTPH=COS(TPH)
	1637	DO I=ITS,MIN(ITE,IDE-1)
	1638
	1639	if (I .eq. ITS) THEN
	1640	TLM=TLM+TDLM*ITS
	1641	else
	1642	TLM=TLM+TDLM
	1643	endif
	1644
	1645	FP=TWOM(CTPH0STPH+STPH0CTPHCOS(TLM))
	1646	F(I,J)=0.5GRID%DTFP
	1647
	1648	ENDDO
	1649	ENDDO
	1650
	1651	! --------------DERIVED VERTICAL GRID CONSTANTS--------------------------
	1652
	1653	EF4T=.5*GRID%DT/CP
	1654	F4Q = -GRID%DT*DTAD
	1655	F4D =-.5GRID%DTDTAD
	1656
	1657	DO L=KDS,KDE-1
	1658	RDETA(L)=1./DETA(L)
	1659	F4Q2(L)=-.25GRID%DTDTAD/DETA(L)
	1660	ENDDO
	1661
	1662	DO J=JTS,min(JTE,JDE-1)
	1663	DO I=ITS,min(ITE,IDE-1)
	1664	DX_NMM(I,J)=DXJ(J)
	1665	WPDAR(I,J)=WPDARJ(J)*HBM2(I,J)
	1666	CPGFU(I,J)=CPGFUJ(J)*VBM2(I,J)
	1667	CURV(I,J)=CURVJ(J)*VBM2(I,J)
	1668	FCP(I,J)=FCPJ(J)*HBM2(I,J)
	1669	FDIV(I,J)=FDIVJ(J)*HBM2(I,J)
	1670	FAD(I,J)=FADJ(J)
	1671	HDACV(I,J)=HDACJ(J)*VBM2(I,J)
	1672	HDAC(I,J)=HDACJ(J)1.25HBM2(I,J)
	1673	ENDDO
	1674	ENDDO
	1675
	1676	DO J=JTS, MIN(JDE-1,JTE)
	1677
	1678	IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN
	1679
	1680	KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have
	1681	DO I=ITS,MIN(IDE-1,ITE)
	1682	IF (I .ge. 2 .and. I .le. KHH) THEN
	1683	HDAC(I,J)=HDAC(I,J)* DFC
	1684	ENDIF
	1685	ENDDO
	1686
	1687	ELSE
	1688
	1689	KHH=2+MOD(J,2)
	1690	DO I=ITS,MIN(IDE-1,ITE)
	1691	IF (I .ge. 2 .and. I .le. KHH) THEN
	1692	HDAC(I,J)=HDAC(I,J)* DFC
	1693	ENDIF
	1694	ENDDO
	1695
	1696	KHH=(IDE-1)-2+MOD(J,2)
	1697
	1698	DO I=ITS,MIN(IDE-1,ITE)
	1699	IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN
	1700	HDAC(I,J)=HDAC(I,J)* DFC
	1701	ENDIF
	1702	ENDDO
	1703	ENDIF
	1704	ENDDO
	1705
	1706	DO J=JTS,min(JTE,JDE-1)
	1707	DO I=ITS,min(ITE,IDE-1)
	1708	DDMPU(I,J)=DDMPUJ(J)*VBM2(I,J)
	1709	DDMPV(I,J)=DDMPVJ(J)*VBM2(I,J)
	1710	HDACV(I,J)=HDACV(I,J)*VBM2(I,J)
	1711	ENDDO
	1712	ENDDO
	1713	! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
	1714
	1715	DO J=JTS,MIN(JDE-1,JTE)
	1716	IF (J.LE.5.OR.J.GE.JDE-1-4) THEN
	1717	KVH=(IDE-1)-1-MOD(J,2)
	1718	DO I=ITS,min(IDE-1,ITE)
	1719	IF (I .ge. 2 .and. I .le. KVH) THEN
	1720	DDMPU(I,J)=DDMPU(I,J)*DDFC
	1721	DDMPV(I,J)=DDMPV(I,J)*DDFC
	1722	HDACV(I,J)=HDACV(I,J)* DFC
	1723	ENDIF
	1724	ENDDO
	1725	ELSE
	1726	KVH=3-MOD(J,2)
	1727	DO I=ITS,min(IDE-1,ITE)
	1728	IF (I .ge. 2 .and. I .le. KVH) THEN
	1729	DDMPU(I,J)=DDMPU(I,J)*DDFC
	1730	DDMPV(I,J)=DDMPV(I,J)*DDFC
	1731	HDACV(I,J)=HDACV(I,J)* DFC
	1732	ENDIF
	1733	ENDDO
	1734	KVH=(IDE-1)-1-MOD(J,2)
	1735	DO I=ITS,min(IDE-1,ITE)
	1736	IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN
	1737	DDMPU(I,J)=DDMPU(I,J)*DDFC
	1738	DDMPV(I,J)=DDMPV(I,J)*DDFC
	1739	HDACV(I,J)=HDACV(I,J)* DFC
	1740	ENDIF
	1741	ENDDO
	1742	ENDIF
	1743	ENDDO
	1744
	1745	write(message,*) 'STC(1)'
	1746	CALL wrf_message(message)
	1747	DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
	1748	write(message,635) (stc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1)
	1749	CALL wrf_message(message)
	1750	ENDDO
	1751
	1752	write(message,*) 'SMC(1)'
	1753	CALL wrf_message(message)
	1754	DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
	1755	write(message,635) (smc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1)
	1756	CALL wrf_message(message)
	1757	ENDDO
	1758
	1759	DO j = jts, MIN(jde-1,jte)
	1760	DO i= ITS, MIN(IDE-1,ITE)
	1761
	1762	if (SM(I,J) .lt. 0.1 .and. SMC(I,1,J) .gt. 0.5 .and. SICE(I,J) .lt. 0.1) then
	1763	write(message,*) 'very moist on land point: ', I,J,SMC(I,1,J)
	1764	CALL wrf_debug(10,message)
	1765	endif
	1766
	1767	enddo
	1768	enddo
	1769
	1770	!!! compute EMT, EM on global domain, and only on task 0.
	1771
	1772	#ifdef DM_PARALLEL
	1773	IF (wrf_dm_on_monitor()) THEN !!!! NECESSARY TO LIMIT THIS TO TASK ZERO?
	1774	#else
	1775	IF (JDS .eq. JTS) THEN !! set unfailable condition for serial job
	1776	#endif
	1777
	1778	ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1))
	1779
	1780	DO J=JDS,JDE-1
	1781	TPH=SB+float(J-1)*DPH
	1782	DXP=ERADDLMCOS(TPH)
	1783	EMJ(J)= GRID%DT/( 4.DXP)DTAD
	1784	EMTJ(J)=GRID%DT/(16.DXP)DTAD
	1785	ENDDO
	1786
	1787	JA=0
	1788	DO 161 J=3,5
	1789	JA=JA+1
	1790	KHLA(JA)=2
	1791	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
	1792	161 EMT(JA)=EMTJ(J)
	1793	DO 162 J=(JDE-1)-4,(JDE-1)-2
	1794	JA=JA+1
	1795	KHLA(JA)=2
	1796	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
	1797	162 EMT(JA)=EMTJ(J)
	1798	DO 163 J=6,(JDE-1)-5
	1799	JA=JA+1
	1800	KHLA(JA)=2
	1801	KHHA(JA)=2+MOD(J,2)
	1802	163 EMT(JA)=EMTJ(J)
	1803	DO 164 J=6,(JDE-1)-5
	1804	JA=JA+1
	1805	KHLA(JA)=(IDE-1)-2
	1806	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
	1807	164 EMT(JA)=EMTJ(J)
	1808
	1809	! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
	1810
	1811	JA=0
	1812	DO 171 J=3,5
	1813	JA=JA+1
	1814	KVLA(JA)=2
	1815	KVHA(JA)=(IDE-1)-1-MOD(J,2)
	1816	171 EM(JA)=EMJ(J)
	1817	DO 172 J=(JDE-1)-4,(JDE-1)-2
	1818	JA=JA+1
	1819	KVLA(JA)=2
	1820	KVHA(JA)=(IDE-1)-1-MOD(J,2)
	1821	172 EM(JA)=EMJ(J)
	1822	DO 173 J=6,(JDE-1)-5
	1823	JA=JA+1
	1824	KVLA(JA)=2
	1825	KVHA(JA)=2+MOD(J+1,2)
	1826	173 EM(JA)=EMJ(J)
	1827	DO 174 J=6,(JDE-1)-5
	1828	JA=JA+1
	1829	KVLA(JA)=(IDE-1)-2
	1830	KVHA(JA)=(IDE-1)-1-MOD(J,2)
	1831	174 EM(JA)=EMJ(J)
	1832
	1833	696 continue
	1834	ENDIF ! wrf_dm_on_monitor/serial job
	1835
	1836	call NMM_SH2O(IMS,IME,JMS,JME,ITS,NNXP,JTS,NNYP,grid%num_soil_layers,ISLTYP, &
	1837	SM,SICE,STC,SMC,SH2O)
	1838
	1839	!! must be a better place to put this, but will eliminate "phantom"
	1840	!! wind points here (no wind point on eastern boundary of odd numbered rows)
	1841
	1842	IF ( abs(IDE-1-ITE) .eq. 1 ) THEN ! along eastern boundary
	1843	write(message,*) 'zero phantom winds'
	1844	CALL wrf_message(message)
	1845	DO K=1,KDE-1
	1846	DO J=JDS,JDE-1,2
	1847	IF (J .ge. JTS .and. J .le. JTE) THEN
	1848	u(IDE-1,J,K)=0.
	1849	v(IDE-1,J,K)=0.
	1850	ENDIF
	1851	ENDDO
	1852	ENDDO
	1853	ENDIF
	1854
	1855	969 continue
	1856
	1857	DO j = jms, jme
	1858	DO i = ims, ime
	1859	fisx=max(fis(i,j),0.)
	1860	Z0(I,J) =SM(I,J)Z0SEA+(1.-SM(I,J)) &
	1861	& (0.Z0MAX+FISx FCM+Z0LAND)
	1862	ENDDO
	1863	ENDDO
	1864
	1865	write(message,*) 'Z0 over memory, leaving module_initialize_real'
	1866	CALL wrf_message(message)
	1867	DO J=JME,JMS,-((JME-JMS)/20+1)
	1868	write(message,635) (Z0(I,J),I=IMS,IME,(IME-IMS)/14+1)
	1869	CALL wrf_message(message)
	1870	ENDDO
	1871
	1872
	1873	endif ! on first_time check
	1874
	1875	write(message,*) 'leaving init_domain_nmm'
	1876	CALL wrf_message( TRIM(message) )
	1877	!
	1878	write(message,*)'STUFF MOVED TO REGISTRY:',grid%IDTAD, &
	1879	& grid%NSOIL,grid%NRADL,grid%NRADS,grid%NPHS,grid%NCNVC,grid%sigma
	1880	CALL wrf_message( TRIM(message) )
	1881	!==================================================================================
	1882
	1883	#define COPY_OUT
	1884	#include <scalar_derefs.inc>
	1885	RETURN
	1886
	1887	END SUBROUTINE init_domain_nmm
	1888
	1889	!------------------------------------------------------
	1890
	1891	SUBROUTINE define_nmm_vertical_coord ( LM, PTSGM, PT, PDTOP,HYBLEVS, &
	1892	SG1,DSG1,SGML1, &
	1893	SG2,DSG2,SGML2,DFL, DFRLG )
	1894
	1895	IMPLICIT NONE
	1896
	1897	! USE module_model_constants
	1898
	1899	!!! certain physical parameters here probably don't need to be defined, as defined
	1900	!!! elsewhere within WRF. Done for initial testing purposes.
	1901
	1902	INTEGER :: LM, LPT2, L
	1903	REAL :: PTSGM, PT, PL, PT2, PDTOP
	1904	REAL :: RGOG, PSIG,PHYB,PHYBM
	1905	REAL, PARAMETER :: Rd = 287.04 ! J deg{-1} kg{-1}
	1906	REAL, PARAMETER :: CP=1004.6,GAMMA=.0065,PRF0=101325.,T0=288.
	1907	REAL, PARAMETER :: g=9.81
	1908
	1909	REAL, DIMENSION(LM) :: DSG,DSG1,DSG2
	1910	REAL, DIMENSION(LM) :: SGML1,SGML2
	1911	REAL, DIMENSION(LM+1) :: SG1,SG2,HYBLEVS,DFL,DFRLG
	1912
	1913	CHARACTER(LEN=255) :: message
	1914
	1915	LPT2=LM+1
	1916
	1917	write(message,*) 'pt= ', pt
	1918	CALL wrf_message(message)
	1919
	1920	DO L=LM+1,1,-1
	1921	pl=HYBLEVS(L)*(101325.-pt)+pt
	1922	if(pl.lt.ptSGm) LPT2=l
	1923	ENDDO
	1924
	1925	IF(LPT2.lt.LM+1) THEN
	1926	pt2=HYBLEVS(LPT2)*(101325.-pt)+pt
	1927	ELSE
	1928	pt2=pt
	1929	ENDIF
	1930
	1931	write(message,) '** Sigma system starts at ',pt2,' Pa, from level ',LPT2
	1932	CALL wrf_message(message)
	1933
	1934	pdtop=pt2-pt
	1935
	1936	write(message,*) 'allocating DSG,DSG1,DSG2 as ', LM
	1937	CALL wrf_debug(10,message)
	1938
	1939	DSG=-99.
	1940
	1941	DO L=1,LM
	1942	DSG(L)=HYBLEVS(L)- HYBLEVS(L+1)
	1943	ENDDO
	1944
	1945	DSG1=0.
	1946	DSG2=0.
	1947
	1948	DO L=LM,1,-1
	1949
	1950	IF(L.ge.LPT2) then
	1951	DSG1(L)=DSG(L)
	1952	ELSE
	1953	DSG2(L)=DSG(L)
	1954	ENDIF
	1955
	1956	ENDDO
	1957
	1958	SGML1=-99.
	1959	SGML2=-99.
	1960
	1961	IF(LPT2.le.LM+1) THEN
	1962
	1963	DO L=LM+1,LPT2,-1
	1964	SG2(L)=0.
	1965	ENDDO
	1966
	1967	DO L=LPT2,2,-1
	1968	SG2(L-1)=SG2(L)+DSG2(L-1)
	1969	ENDDO
	1970
	1971	DO L=LPT2-1,1,-1
	1972	SG2(L)=SG2(L)/SG2(1)
	1973	ENDDO
	1974	SG2(1)=1.
	1975
	1976	DO L=LPT2-1,1,-1
	1977	DSG2(L)=SG2(L)-SG2(L+1)
	1978	SGML2(l)=(SG2(l)+SG2(l+1))*0.5
	1979	ENDDO
	1980
	1981	ENDIF
	1982
	1983	DO L=LM,LPT2,-1
	1984	DSG2(L)=0.
	1985	SGML2(L)=0.
	1986	ENDDO
	1987
	1988	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1989
	1990	SG1(LM+1)=0.
	1991
	1992	DO L=LM+1,LPT2,-1
	1993	SG1(L-1)=SG1(L)+DSG1(L-1)
	1994	ENDDO
	1995
	1996	DO L=LM,LPT2,-1
	1997	SG1(L)=SG1(L)/SG1(LPT2-1)
	1998	ENDDO
	1999
	2000	SG1(LPT2-1)=1.
	2001
	2002	do l=LPT2-2,1,-1
	2003	SG1(l)=1.
	2004	enddo
	2005
	2006
	2007	DO L=LM,LPT2,-1
	2008	DSG1(L)=SG1(L)-SG1(L+1)
	2009	SGML1(L)=(SG1(L)+SG1(L+1))*0.5
	2010	ENDDO
	2011
	2012	DO L=LPT2-1,1,-1
	2013	DSG1(L)=0.
	2014	SGML1(L)=1.
	2015	ENDDO
	2016
	2017	1000 format('l,hyblevs,psig,SG1,SG2,phyb,phybm')
	2018	1100 format(' ',i4,f7.4,f10.2,2f7.4,2f10.2)
	2019
	2020	write(message,1000)
	2021	CALL wrf_debug(100,message)
	2022
	2023	do l=1,LM+1
	2024	psig=HYBLEVS(L)*(101325.-pt)+pt
	2025	phyb=SG1(l)pdtop+SG2(l)(101325.-pdtop-pt)+pt
	2026	if(l.lt.LM+1) then
	2027	phybm=SGML1(l)pdtop+SGML2(l)(101325.-pdtop-pt)+pt
	2028	else
	2029	phybm=-99.
	2030	endif
	2031
	2032	write(message,1100) l,HYBLEVS(L),psig &
	2033	,SG1(l),SG2(l),phyb,phybm
	2034	CALL wrf_debug(100,message)
	2035	enddo
	2036
	2037
	2038	632 format(f9.6)
	2039
	2040	write(message,*) 'SG1'
	2041	CALL wrf_debug(100,message)
	2042	do L=LM+1,1,-1
	2043	write(message,632) SG1(L)
	2044	CALL wrf_debug(100,message)
	2045	enddo
	2046
	2047	write(message,*) 'SG2'
	2048	CALL wrf_debug(100,message)
	2049	do L=LM+1,1,-1
	2050	write(message,632) SG2(L)
	2051	CALL wrf_debug(100,message)
	2052	enddo
	2053
	2054	write(message,*) 'DSG1'
	2055	CALL wrf_debug(100,message)
	2056	do L=LM,1,-1
	2057	write(message,632) DSG1(L)
	2058	CALL wrf_debug(100,message)
	2059	enddo
	2060
	2061	write(message,*) 'DSG2'
	2062	CALL wrf_debug(100,message)
	2063	do L=LM,1,-1
	2064	write(message,632) DSG2(L)
	2065	CALL wrf_debug(100,message)
	2066	enddo
	2067
	2068	write(message,*) 'SGML1'
	2069	CALL wrf_debug(100,message)
	2070	do L=LM,1,-1
	2071	write(message,632) SGML1(L)
	2072	CALL wrf_debug(100,message)
	2073	enddo
	2074
	2075	write(message,*) 'SGML2'
	2076	CALL wrf_debug(100,message)
	2077	do L=LM,1,-1
	2078	write(message,632) SGML2(L)
	2079	CALL wrf_debug(100,message)
	2080	enddo
	2081
	2082	rgog=(rd*gamma)/g
	2083	DO L=1,LM+1
	2084	DFL(L)=gT0(1.-((pt+SG1(L)pdtop+SG2(L)(101325.-pt2)) &
	2085	/101325.)**rgog)/gamma
	2086	DFRLG(L)=DFL(L)/g
	2087	write(message,*) 'L, DFL(L): ', L, DFL(L)
	2088	CALL wrf_debug(100,message)
	2089	ENDDO
	2090
	2091	END SUBROUTINE define_nmm_vertical_coord
	2092
	2093	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2094	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2095	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2096
	2097	SUBROUTINE compute_nmm_surfacep ( TERRAIN_HGT_T, Z3D_IN, PRESS3D_IN, T3D_IN &
	2098	&, psfc_out,generic &
	2099	&, IDS,IDE,JDS,JDE,KDS,KDE &
	2100	&, IMS,IME,JMS,JME,KMS,KME &
	2101	&, ITS,ITE,JTS,JTE,KTS,KTE )
	2102
	2103
	2104	IMPLICIT NONE
	2105
	2106	real, allocatable:: dum2d(:,:),DUM2DB(:,:)
	2107
	2108	integer :: IDS,IDE,JDS,JDE,KDS,KDE
	2109	integer :: IMS,IME,JMS,JME,KMS,KME
	2110	integer :: ITS,ITE,JTS,JTE,KTS,KTE,Ilook,Jlook
	2111	integer :: I,J,II,generic,L,KINSERT,K,bot_lev,LL
	2112	integer :: IHE(JMS:JME),IHW(JMS:JME), loopinc,iloopinc
	2113
	2114	real :: TERRAIN_HGT_T(IMS:IME,JMS:JME)
	2115	real :: Z3D_IN(IMS:IME,JMS:JME,generic)
	2116	real :: T3D_IN(IMS:IME,JMS:JME,generic)
	2117	real :: PRESS3D_IN(IMS:IME,JMS:JME,generic)
	2118	real :: PSFC_IN(IMS:IME,JMS:JME),TOPO_IN(IMS:IME,JMS:JME)
	2119	real :: psfc_out(IMS:IME,JMS:JME),rincr(IMS:IME,JMS:JME)
	2120	real :: dif1,dif2,dif3,dif4,dlnpdz,BOT_INPUT_HGT,BOT_INPUT_PRESS,dpdz,rhs
	2121	real :: zin(generic),pin(generic)
	2122
	2123	character (len=255) :: message
	2124
	2125	logical :: DEFINED_PSFC(IMS:IME,JMS:JME), DEFINED_PSFCB(IMS:IME,JMS:JME)
	2126
	2127	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2128
	2129	Ilook=25
	2130	Jlook=25
	2131
	2132	DO j = JMS, JME
	2133	IHE(J)=MOD(J+1,2)
	2134	IHW(J)=IHE(J)-1
	2135	ENDDO
	2136
	2137	DO J=JMS,JME
	2138	DO I=IMS,IME
	2139	DEFINED_PSFC(I,J)=.FALSE.
	2140	DEFINED_PSFCB(I,J)=.FALSE.
	2141	IF (PRESS3D_IN(I,J,1) .ne. 200100.) THEN
	2142	PSFC_IN(I,J)=PRESS3D_IN(I,J,1)
	2143	TOPO_IN(I,J)=Z3D_IN(I,J,1)
	2144	ELSE
	2145	PSFC_IN(I,J)=PRESS3D_IN(I,J,2)
	2146	TOPO_IN(I,J)=Z3D_IN(I,J,2)
	2147	ENDIF
	2148	ENDDO
	2149	ENDDO
	2150
	2151	! input surface pressure smoothing over the ocean - still needed for NAM?
	2152
	2153	II_loop: do II=1,8
	2154
	2155	CYCLE II_loop
	2156
	2157	do J=JTS+1,min(JTE,JDE-1)-1
	2158	do I=ITS+1,min(ITE,IDE-1)-1
	2159	rincr(I,J)=0.
	2160
	2161	if (PSFC_IN(I,J) .gt. 100000. .and. &
	2162	PSFC_IN(I+IHE(J),J+1) .gt. 100000. .and. &
	2163	PSFC_IN(I+IHE(J),J-1) .gt. 100000. .and. &
	2164	PSFC_IN(I+IHW(J),J+1) .gt. 100000. .and. &
	2165	PSFC_IN(I+IHW(J),J-1) .gt. 100000. ) then
	2166
	2167	dif1=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J+1))
	2168	dif2=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J-1))
	2169	dif3=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J+1))
	2170	dif4=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J-1))
	2171
	2172	if (max(dif1,dif2,dif3,dif4) .lt. 200. .and. TOPO_IN(I,J).le. 0.5 .and. &
	2173	TOPO_IN(I+IHE(J),J+1) .le. 0.5 .and. &
	2174	TOPO_IN(I+IHW(J),J+1) .le. 0.5 .and. &
	2175	TOPO_IN(I+IHE(J),J-1) .le. 0.5 .and. &
	2176	TOPO_IN(I+IHW(J),J-1) .lt. 0.5) then
	2177
	2178	rincr(I,J)=0.125( 4.PSFC_IN(I,J)+ &
	2179	PSFC_IN(I+IHE(J),J+1)+PSFC_IN(I+IHE(J),J-1)+ &
	2180	PSFC_IN(I+IHW(J),J+1)+PSFC_IN(I+IHW(J),J-1) ) &
	2181	- PSFC_IN(I,J)
	2182
	2183	! if (rincr(I,J) .ne. 0 .and. abs(rincr(I,J)) .gt. 20.) then
	2184	! write(message,*) 'II, I,J,rincr: ', II, I,J,rincr(I,J)
	2185	! CALL wrf_message(message)
	2186	! endif
	2187
	2188	endif
	2189	endif
	2190
	2191	ENDDO
	2192	ENDDO
	2193
	2194	DO J=JTS+1,min(JTE,JDE-1)-1
	2195	DO I=ITS+1,min(ITE,IDE-1)-1
	2196	PSFC_IN(I,J)=PSFC_IN(I,J) + rincr(I,J)
	2197	ENDDO
	2198	ENDDO
	2199
	2200	! write(message,*) ' -------------------------------------------------- '
	2201	! CALL wrf_message(message)
	2202
	2203	end do II_loop
	2204
	2205	ALLOCATE(DUM2D(IMS:IME,JMS:JME))
	2206
	2207	DO J=JMS,JME
	2208	DO I=IMS,IME
	2209	DUM2D(I,J)=-9.
	2210	END DO
	2211	END DO
	2212
	2213	DO J=JTS,min(JTE,JDE-1)
	2214	I_loop: DO I=ITS,min(ITE,IDE-1)
	2215
	2216	IF (PSFC_IN(I,J) .lt. 0.1) THEN
	2217	write(message,*) 'QUITTING BECAUSE I,J, PSFC_IN: ', I,J,PSFC_IN(I,J)
	2218	call wrf_error_fatal(message)
	2219	ENDIF
	2220
	2221	BOT_INPUT_PRESS=PSFC_IN(I,J)
	2222	BOT_INPUT_HGT=TOPO_IN(I,J)
	2223
	2224	IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
	2225
	2226	write(message,*) ' TERRAIN_HGT_T: ', I,J, TERRAIN_HGT_T(I,J)
	2227	CALL wrf_message(message)
	2228	write(message,*) ' PSFC_IN, TOPO_IN: ', &
	2229	I, J, PSFC_IN(I,J),TOPO_IN(I,J)
	2230	CALL wrf_message(message)
	2231
	2232	DO L=1,generic
	2233	write(message,*) ' L,PRESS3D_IN, Z3D_IN: ', &
	2234	I,J,L, PRESS3D_IN(I,J,L),Z3D_IN(I,J,L)
	2235	CALL wrf_debug(10,message)
	2236	END DO
	2237	ENDIF
	2238
	2239	DO L=2,generic-1
	2240
	2241	IF ( PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. &
	2242	Z3D_IN(I,J,L) .lt. TERRAIN_HGT_T(I,J) .AND. &
	2243	Z3D_IN(I,J,L+1) .gt. TERRAIN_HGT_T(I,J) ) THEN
	2244
	2245	BOT_INPUT_PRESS=PRESS3D_IN(i,j,L)
	2246	BOT_INPUT_HGT=Z3D_IN(I,J,L)
	2247
	2248	! IF (I .eq. Ilook .and. J .eq. Jlook) THEN
	2249	! write(message,*) 'BOT_INPUT_PRESS, BOT_INPUT_HGT NOW : ', &
	2250	! Ilook,Jlook, BOT_INPUT_PRESS, BOT_INPUT_HGT
	2251	! CALL wrf_message(message)
	2252	! ENDIF
	2253
	2254	ENDIF
	2255	END DO
	2256
	2257	!!!!!!!!!!!!!!!!!!!!!! START HYDRO CHECK
	2258
	2259	IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. &
	2260	(PSFC_IN(I,J) .gt. PRESS3D_IN(i,j,2) .OR. &
	2261	TOPO_IN(I,J) .lt. Z3D_IN(I,J,2)) ) THEN ! extrapolate downward
	2262
	2263	IF (J .eq. JTS .AND. I .eq. ITS) THEN
	2264	write(message,*) 'hydro check - should only be for isobaric input'
	2265	CALL wrf_message(message)
	2266	ENDIF
	2267
	2268	IF (Z3D_IN(I,J,2) .ne. TOPO_IN(I,J)) THEN
	2269	dpdz=(PRESS3D_IN(i,j,2)-PSFC_IN(I,J))/(Z3D_IN(I,J,2)-TOPO_IN(I,J))
	2270	rhs=-9.81((PRESS3D_IN(i,j,2)+ PSFC_IN(I,J))/2.)/(287.04 T3D_IN(I,J,2))
	2271
	2272	IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN
	2273	IF (dpdz .lt. 1.05rhs .OR. dpdz .gt. 0.95rhs) THEN
	2274	write(message,*) 'I,J,P(2),Psfc,Z(2),Zsfc: ', &
	2275	I,J,PRESS3D_IN(i,j,2),PSFC_IN(I,J),Z3D_IN(I,J,2),TOPO_IN(I,J)
	2276	IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) CALL wrf_debug(50,message)
	2277	CYCLE I_loop
	2278	ENDIF
	2279
	2280	ENDIF
	2281
	2282	ELSE ! z(2) equals TOPO_IN
	2283
	2284	IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN
	2285	! write(message,*) 'all equal at I,J: ', I,J
	2286	! CALL wrf_message(message)
	2287	ELSE
	2288	! write(message,*) 'heights equal, pressures not: ', &
	2289	! PRESS3D_IN(i,j,2), PSFC_IN(I,J)
	2290	! CALL wrf_message(message)
	2291	CYCLE I_loop
	2292	ENDIF
	2293
	2294	ENDIF
	2295
	2296	IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN
	2297	IF (PRESS3D_IN(i,j,2) .lt. PSFC_IN(I,J) .and. &
	2298	Z3D_IN(I,J,2) .lt. TOPO_IN(I,J)) THEN
	2299	! write(message,*) 'surface data mismatch(a) at I,J: ', I,J
	2300	! CALL wrf_message(message)
	2301	CYCLE I_loop
	2302	ELSEIF (PRESS3D_IN(i,j,2) .gt. PSFC_IN(I,J) .AND. &
	2303	Z3D_IN(I,J,2) .gt. TOPO_IN(I,J)) THEN
	2304	! write(message,*) 'surface data mismatch(b) at I,J: ', I,J
	2305	! CALL wrf_message(message)
	2306	CYCLE I_loop
	2307	ENDIF
	2308	ENDIF
	2309	ENDIF
	2310
	2311	!!!!!!! loop over a few more levels
	2312
	2313	DO L=3,6
	2314	IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. &
	2315	(((PSFC_IN(I,J)-PRESS3D_IN(i,j,L)) .lt. 400.) .OR. &
	2316	TOPO_IN(I,J) .lt. Z3D_IN(I,J,L))) then
	2317
	2318	IF (Z3D_IN(I,J,L) .ne. TOPO_IN(I,J)) THEN
	2319	dpdz=(PRESS3D_IN(i,j,L)-PSFC_IN(I,J))/ &
	2320	(Z3D_IN(I,J,L)-TOPO_IN(I,J))
	2321	rhs=-9.81*((PRESS3D_IN(i,j,L)+ PSFC_IN(I,J))/2.)/ &
	2322	(287.04*T3D_IN(I,J,L))
	2323	IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN
	2324	IF (dpdz .lt. 1.05rhs .or. dpdz .gt. 0.95rhs) THEN
	2325	write(message,*) 'I,J,L,Piso,Psfc,Ziso,Zsfc: ', &
	2326	I,J,L,PRESS3D_IN(i,j,L),PSFC_IN(I,J),&
	2327	Z3D_IN(I,J,L),TOPO_IN(I,J)
	2328	IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) &
	2329	CALL wrf_debug(50,message)
	2330	CYCLE I_loop
	2331	ENDIF
	2332	ENDIF
	2333	ELSE
	2334	IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN
	2335	! write(message,*) 'all equal at I,J: ', I,J
	2336	! CALL wrf_message(message)
	2337	ELSE
	2338	CYCLE I_loop
	2339	ENDIF
	2340	ENDIF
	2341	ENDIF
	2342
	2343	IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN
	2344	IF (PRESS3D_IN(i,j,L) .lt. PSFC_IN(I,J) .AND. &
	2345	Z3D_IN(I,J,L) .lt. TOPO_IN(I,J)) THEN
	2346	CYCLE I_loop
	2347	ELSEIF (PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. &
	2348	Z3D_IN(I,J,L) .gt. TOPO_IN(I,J)) THEN
	2349	CYCLE I_loop
	2350	ENDIF
	2351	ENDIF
	2352	END DO
	2353	!!!!!!!!!!!!!!!!!!!!!! END HYDRO CHECK
	2354
	2355	IF (TERRAIN_HGT_T(I,J) .eq. BOT_INPUT_HGT ) THEN
	2356	dum2d(I,J)=BOT_INPUT_PRESS
	2357
	2358	IF (BOT_INPUT_HGT .ne. 0. .and. (BOT_INPUT_HGT-INT(BOT_INPUT_HGT) .ne. 0.) ) THEN
	2359	write(message,*) 'with BOT_INPUT_HGT: ', BOT_INPUT_HGT, &
	2360	'set dum2d to bot_input_pres: ', I,J,dum2d(I,J)
	2361	CALL wrf_message(message)
	2362	ENDIF
	2363
	2364	IF (dum2d(I,J) .lt. 50000. .OR. dum2d(I,J) .gt. 109000.) THEN
	2365	write(message,*) 'bad dum2d(a): ', I,J,DUM2D(I,J)
	2366	CALL wrf_message(message)
	2367	ENDIF
	2368
	2369	ELSEIF (TERRAIN_HGT_T(I,J) .lt. BOT_INPUT_HGT ) THEN
	2370
	2371	! target is below lowest possible input...extrapolate
	2372
	2373	IF ( BOT_INPUT_PRESS-PRESS3D_IN(I,J,2) .gt. 500. ) THEN
	2374	dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,2)) ) / &
	2375	(BOT_INPUT_HGT-Z3D_IN(i,j,2))
	2376	IF (I .eq. Ilook .and. J .eq. Jlook) THEN
	2377	write(message,*) 'I,J,dlnpdz(a): ', I,J,dlnpdz
	2378	CALL wrf_message(message)
	2379	ENDIF
	2380
	2381	ELSE
	2382
	2383	!! thin layer and/or just have lowest level - difference with 3rd level data
	2384	IF ( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,3)) .gt. 290. ) THEN
	2385
	2386	dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,3)) ) / &
	2387	(BOT_INPUT_HGT-Z3D_IN(i,j,3))
	2388
	2389	IF (I .eq. Ilook .and. J .eq. Jlook) then
	2390	write(message,*) 'p diff: ', BOT_INPUT_PRESS, PRESS3D_IN(i,j,3)
	2391	CALL wrf_message(message)
	2392	write(message,*) 'z diff: ', BOT_INPUT_HGT, Z3D_IN(i,j,3)
	2393	CALL wrf_message(message)
	2394	ENDIF
	2395
	2396	ELSE
	2397
	2398	!! Loop up to level 7 looking for a sufficiently thick layer
	2399
	2400	FIND_THICK: DO LL=4,7
	2401	IF( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,LL)) .gt. 290.) THEN
	2402	dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,LL)) ) / &
	2403	(BOT_INPUT_HGT-Z3D_IN(i,j,LL))
	2404	EXIT FIND_THICK
	2405	ENDIF
	2406	END DO FIND_THICK
	2407
	2408	ENDIF
	2409
	2410	ENDIF
	2411
	2412	dum2d(I,J)= exp(log(BOT_INPUT_PRESS) + dlnpdz * &
	2413	(TERRAIN_HGT_T(I,J) - BOT_INPUT_HGT) )
	2414
	2415	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
	2416	write(message,*) 'bad dum2d(b): ', I,J,DUM2D(I,J)
	2417	CALL wrf_message(message)
	2418	write(message,*) 'BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T, BOT_INPUT_HGT: ', &
	2419	BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T(I,J), BOT_INPUT_HGT
	2420	CALL wrf_message(message)
	2421	write(message,*) 'Z3D_IN: ', Z3D_IN(I,J,1:10)
	2422	CALL wrf_message(message)
	2423	write(message,*) 'PRESS3D_IN: ', PRESS3D_IN(I,J,1:10)
	2424	CALL wrf_message(message)
	2425	ENDIF
	2426
	2427	ELSE ! target level bounded by input levels
	2428
	2429	DO L=2,generic-1
	2430	IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. &
	2431	TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN
	2432	dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / &
	2433	(Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1))
	2434	dum2d(I,J)= log(PRESS3D_IN(i,j,l)) + &
	2435	dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) )
	2436	dum2d(i,j)=exp(dum2d(i,j))
	2437	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
	2438	write(message,*) 'bad dum2d(c): ', I,J,DUM2D(I,J)
	2439	CALL wrf_message(message)
	2440	ENDIF
	2441	ENDIF
	2442	ENDDO
	2443
	2444	!!! account for situation where BOT_INPUT_HGT < TERRAIN_HGT_T < Z3D_IN(:,2,:)
	2445	IF (dum2d(I,J) .eq. -9 .AND. BOT_INPUT_HGT .lt. TERRAIN_HGT_T(I,J) &
	2446	.AND. TERRAIN_HGT_T(I,J) .lt. Z3D_IN(I,J,2)) then
	2447
	2448	IF (mod(I,50) .eq. 0 .AND. mod(J,50) .eq. 0) THEN
	2449	write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', &
	2450	I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J)
	2451	CALL wrf_message(message)
	2452	ENDIF
	2453
	2454	dlnpdz= (log(PSFC_IN(i,j))-log(PRESS3D_IN(i,j,2)) ) / &
	2455	(TOPO_IN(i,j)-Z3D_IN(i,j,2))
	2456	dum2d(I,J)= log(PSFC_IN(i,j)) + &
	2457	dlnpdz * (TERRAIN_HGT_T(I,J) - TOPO_IN(i,j) )
	2458	dum2d(i,j)= exp(dum2d(i,j))
	2459	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
	2460	write(message,*) 'bad dum2d(d): ', I,J,DUM2D(I,J)
	2461	CALL wrf_message(message)
	2462	ENDIF
	2463	ENDIF
	2464
	2465	IF (dum2d(I,J) .eq. -9.) THEN
	2466	write(message,*) 'must have flukey situation in new ', I,J
	2467	CALL wrf_message(message)
	2468	write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', &
	2469	I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J)
	2470	CALL wrf_message(message)
	2471
	2472	DO L=1,generic-1
	2473	IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN
	2474	! problematic with HGT_M substitution for "input" surface height?
	2475	dum2d(i,j)=PRESS3D_IN(I,J,L)
	2476	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
	2477	write(message,*) 'bad dum2d(e): ', I,J,DUM2D(I,J)
	2478	CALL wrf_message(message)
	2479	ENDIF
	2480	ENDIF
	2481	ENDDO
	2482
	2483	IF ( TERRAIN_HGT_T(I,J) .eq. TOPO_IN(I,J)) THEN
	2484	dum2d(I,J)=PSFC_IN(I,J)
	2485	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
	2486	write(message,*) 'bad dum2d(f): ', I,J,DUM2D(I,J)
	2487	CALL wrf_message(message)
	2488	ENDIF
	2489	write(message,*) 'matched input topo, psfc: ', I,J,TOPO_IN(I,J),PSFC_IN(I,J)
	2490	CALL wrf_message(message)
	2491	ENDIF
	2492
	2493	IF (dum2d(I,J) .eq. -9.) THEN
	2494	CALL wrf_error_fatal("quitting due to undefined surface pressure")
	2495	ENDIF
	2496	ENDIF
	2497
	2498	DEFINED_PSFC(I,J)=.TRUE.
	2499
	2500	IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
	2501	write(message,*) 'newstyle psfc: ', I,J,dum2d(I,J)
	2502	CALL wrf_message(message)
	2503	ENDIF
	2504
	2505	ENDIF
	2506
	2507	ENDDO I_loop
	2508	ENDDO
	2509
	2510	write(message,*) 'psfc points (new style)'
	2511	CALL wrf_message(message)
	2512	loopinc=max( (JTE-JTS)/20,1)
	2513	iloopinc=max( (ITE-ITS)/10,1)
	2514
	2515	DO J=min(JTE,JDE-1),JTS,-loopinc
	2516	write(message,633) (dum2d(I,J)/100.,I=ITS,min(ITE,IDE-1),iloopinc)
	2517	END DO
	2518
	2519	633 format(35(f5.0,1x))
	2520
	2521	write(message,*) 'PSFC extremes (new style)'
	2522	CALL wrf_message(message)
	2523	write(message,*) minval(dum2d,MASK=DEFINED_PSFC),maxval(dum2d,MASK=DEFINED_PSFC)
	2524	CALL wrf_message(message)
	2525
	2526	! IF (minval(dum2d,MASK=DEFINED_PSFC) .lt. 50000. .or. maxval(dum2d,MASK=DEFINED_PSFC) .gt. 108000.) THEN
	2527
	2528	DO J=JTS,min(JTE,JDE-1)
	2529	DO I=ITS,min(ITE,IDE-1)
	2530
	2531	IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .lt. 50000. ) THEN
	2532	IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
	2533	WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
	2534	CALL wrf_debug(2,message)
	2535	ELSE
	2536	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
	2537	ENDIF
	2538	ENDIF
	2539
	2540	IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .gt. 108000. ) THEN
	2541	IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
	2542	WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
	2543	CALL wrf_debug(2,message)
	2544	ELSE
	2545	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
	2546	ENDIF
	2547	ENDIF
	2548
	2549	END DO
	2550	END DO
	2551
	2552
	2553
	2554	!! "traditional" isobaric only approach ------------------------------------------------
	2555
	2556	ALLOCATE (DUM2DB(IMS:IME,JMS:JME))
	2557	DO J=JMS,JME
	2558	DO I=IMS,IME
	2559	DUM2DB(I,J)=-9.
	2560	END DO
	2561	END DO
	2562
	2563	DO J=JTS,min(JTE,JDE-1)
	2564	DO I=ITS,min(ITE,IDE-1)
	2565
	2566	IF (TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,2)) THEN ! targ below lowest
	2567
	2568	IF ( abs(PRESS3D_IN(i,j,2)-PRESS3D_IN(i,j,3)) .gt. 290.) THEN
	2569	dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,3)) ) / &
	2570	(Z3D_IN(i,j,2)-Z3D_IN(i,j,3))
	2571	ELSE
	2572	dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,4)) ) / &
	2573	(Z3D_IN(i,j,2)-Z3D_IN(i,j,4))
	2574	ENDIF
	2575
	2576	DUM2DB(I,J)= exp( log(PRESS3D_IN(i,j,2)) + dlnpdz * &
	2577	(TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,2)) )
	2578
	2579	IF (I .eq. Ilook .and. J .eq. Jlook) THEN
	2580	write(message,*) 'I,K, trad: dlnpdz, press_in(2), terrain_t, Z3D_IN(2): ', I,J,dlnpdz, &
	2581	PRESS3D_IN(i,j,2), TERRAIN_HGT_T(I,J), Z3D_IN(i,j,2)
	2582	CALL wrf_message(message)
	2583	ENDIF
	2584
	2585	DEFINED_PSFCB(i,j)=.true.
	2586
	2587	ELSEIF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,2)) THEN ! target level bounded by input levels
	2588
	2589	DO L=2,generic-1
	2590	IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. &
	2591	TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN
	2592
	2593	dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / &
	2594	(Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1))
	2595
	2596	DUM2DB(I,J)= log(PRESS3D_IN(i,j,l)) + &
	2597	dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) )
	2598	DUM2DB(i,j)=exp(DUM2DB(i,j))
	2599
	2600	DEFINED_PSFCB(i,j)=.true.
	2601
	2602	IF (DUM2DB(I,J) .lt. 13000.) THEN
	2603	write(message,*) 'I,J,L,terrain,Z3d(L),z3d(L+1),p3d(L),p3d(l+1): ', I,J,L, &
	2604	TERRAIN_HGT_T(I,J),Z3D_IN(I,J,L),Z3D_IN(I,J,L+1),PRESS3D_IN(I,J,L), &
	2605	PRESS3D_IN(I,J,L+1)
	2606	CALL wrf_error_fatal(message)
	2607	ENDIF
	2608	ENDIF
	2609	ENDDO
	2610
	2611	ELSEIF (TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,2)) THEN
	2612	DUM2DB(i,j)=PRESS3D_IN(I,J,2)
	2613	DEFINED_PSFCB(i,j)=.true.
	2614	ENDIF
	2615
	2616	IF (DUM2DB(I,J) .eq. -9.) THEN
	2617	write(message,*) 'must have flukey situation in trad ', I,J
	2618	CALL wrf_message(message)
	2619	DO L=1,generic-1
	2620	IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN
	2621	DUM2DB(i,j)=PRESS3D_IN(I,J,L)
	2622	DEFINED_PSFCB(i,j)=.true.
	2623	ENDIF
	2624	ENDDO
	2625	ENDIF
	2626
	2627	IF (DUM2DB(I,J) .eq. -9.) THEN
	2628	write(message,*) 'HOPELESS PSFC, I QUIT'
	2629	CALL wrf_error_fatal(message)
	2630	ENDIF
	2631
	2632	if (I .eq. Ilook .and. J .eq. Jlook) THEN
	2633	write(message,*) ' traditional psfc: ', I,J,DUM2DB(I,J)
	2634	CALL wrf_message(message)
	2635	ENDIF
	2636
	2637	ENDDO
	2638	ENDDO
	2639
	2640	write(message,*) 'psfc points (traditional)'
	2641	CALL wrf_message(message)
	2642	DO J=min(JTE,JDE-1),JTS,-loopinc
	2643	write(message,633) (DUM2DB(I,J)/100.,I=its,min(ite,IDE-1),iloopinc)
	2644	CALL wrf_message(message)
	2645	ENDDO
	2646
	2647	write(message,*) 'PSFC extremes (traditional)'
	2648	CALL wrf_message(message)
	2649	write(message,*) minval(DUM2DB,MASK=DEFINED_PSFCB),maxval(DUM2DB,MASK=DEFINED_PSFCB)
	2650	CALL wrf_message(message)
	2651
	2652	DO J=JTS,min(JTE,JDE-1)
	2653	DO I=ITS,min(ITE,IDE-1)
	2654
	2655	IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .lt. 50000. ) THEN
	2656	IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
	2657	WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
	2658	CALL wrf_debug(2,message)
	2659	ELSE
	2660	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
	2661	ENDIF
	2662	ENDIF
	2663
	2664	IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .gt. 108000. ) THEN
	2665	IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
	2666	WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
	2667	CALL wrf_debug(2,message)
	2668	ELSE
	2669	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
	2670	ENDIF
	2671	ENDIF
	2672
	2673	! IF (DEFINED_PSFCB(I,J) .AND. ( dum2db(I,J) .lt. 50000. .OR. dum2db(I,J) .gt. 108000. )) THEN
	2674	! IF (TERRAIN_HGT_T(I,J) .gt. -10. .and. TERRAIN_HGT_T(I,J) .lt. 5000.) THEN
	2675	! write(0,*) 'I,J, terrain_hgt_t, dum2db: ', I,J, terrain_hgt_t(I,J),dum2db(I,J)
	2676	! CALL wrf_error_fatal("quit due to unrealistic surface pressure")
	2677	! ELSE
	2678	! WRITE(message,*) 'surface pressure allowed because surface height is extreme value of: ', TERRAIN_HGT_T(I,J)
	2679	! CALL wrf_debug(2,message)
	2680	! ENDIF
	2681	! ENDIF
	2682
	2683	ENDDO
	2684	ENDDO
	2685
	2686	!!!!! end traditional
	2687
	2688	DO J=JTS,min(JTE,JDE-1)
	2689	DO I=ITS,min(ITE,IDE-1)
	2690	IF (DEFINED_PSFCB(I,J) .and. DEFINED_PSFC(I,J)) THEN
	2691
	2692	IF ( abs(dum2d(I,J)-DUM2DB(I,J)) .gt. 400.) THEN
	2693	write(message,*) 'BIG DIFF I,J, dum2d, DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J)
	2694	CALL wrf_message(message)
	2695	ENDIF
	2696
	2697	!! do we have enough confidence in new style to give it more than 50% weight?
	2698	psfc_out(I,J)=0.5*(dum2d(I,J)+DUM2DB(I,J))
	2699
	2700	ELSEIF (DEFINED_PSFC(I,J)) THEN
	2701	psfc_out(I,J)=dum2d(I,J)
	2702	ELSEIF (DEFINED_PSFCB(I,J)) THEN
	2703	psfc_out(I,J)=DUM2DB(I,J)
	2704	ELSE
	2705	write(message,*) 'I,J,dum2d,DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J)
	2706	CALL wrf_message(message)
	2707	write(message,*) 'I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J): ', I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J)
	2708	CALL wrf_message(message)
	2709	call wrf_error_fatal("psfc_out completely undefined")
	2710	ENDIF
	2711
	2712	IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
	2713	write(message,*) ' combined psfc: ', I,J,psfc_out(I,J)
	2714	CALL wrf_message(message)
	2715	ENDIF
	2716
	2717	IF (PSFC_OUT(I,J) .lt. 50000. ) THEN
	2718	IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
	2719	WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
	2720	CALL wrf_debug(2,message)
	2721	ELSE
	2722	write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J)
	2723	CALL wrf_debug(2,message)
	2724	write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J)
	2725	CALL wrf_debug(2,message)
	2726	write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J)
	2727	CALL wrf_debug(2,message)
	2728	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
	2729	ENDIF
	2730	ENDIF
	2731
	2732	IF (PSFC_OUT(I,J) .gt. 108000. ) THEN
	2733	IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
	2734	WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
	2735	CALL wrf_debug(2,message)
	2736	ELSE
	2737	write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J)
	2738	CALL wrf_debug(2,message)
	2739	write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J)
	2740	CALL wrf_debug(2,message)
	2741	write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J)
	2742	CALL wrf_debug(2,message)
	2743	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
	2744	ENDIF
	2745	ENDIF
	2746
	2747	ENDDO
	2748	ENDDO
	2749
	2750	deallocate(dum2d,dum2db)
	2751
	2752	END SUBROUTINE compute_nmm_surfacep
	2753
	2754	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2755	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2756
	2757	SUBROUTINE compute_3d_pressure(psfc_out,SGML1,SGML2,pdtop,pt &
	2758	&, pd,p3d_out &
	2759	&, IDS,IDE,JDS,JDE,KDS,KDE &
	2760	&, IMS,IME,JMS,JME,KMS,KME &
	2761	&, ITS,ITE,JTS,JTE,KTS,KTE )
	2762
	2763
	2764	INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
	2765	INTEGER :: IMS,IME,JMS,JME,KMS,KME
	2766	INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE
	2767
	2768	REAL, INTENT(IN) :: psfc_out(IMS:IME,JMS:JME)
	2769	REAL, INTENT(IN) :: SGML1(KDE),SGML2(KDE),pdtop,pt
	2770
	2771	REAL, INTENT(OUT):: p3d_out(IMS:IME,JMS:JME,KDS:KDE-1)
	2772	REAL, INTENT(OUT):: PD(IMS:IME,JMS:JME)
	2773
	2774	CHARACTER (len=255) :: message
	2775
	2776	! write(message,*) 'pdtop, pt, psfc_out(1,1): ', pdtop, pt, psfc_out(1,1)
	2777	! CALL wrf_message(message)
	2778
	2779	DO J=JTS,min(JTE,JDE-1)
	2780	DO I=ITS,min(ITE,IDE-1)
	2781	PD(I,J)=psfc_out(I,J)-PDTOP-PT
	2782	ENDDO
	2783	ENDDO
	2784
	2785	DO J=JTS,min(JTE,JDE-1)
	2786	DO K=KDS,KDE-1
	2787	DO I=ITS,min(ITE,IDE-1)
	2788	p3d_out(I,J,K)=PD(I,J)SGML2(K)+PDTOPSGML1(K)+PT
	2789
	2790	IF (p3d_out(I,J,K) .ge. psfc_out(I,J) .or. p3d_out(I,J,K) .le. pt) THEN
	2791	write(message,*) 'I,K,J,p3d_out: ', I,K,J,p3d_out(I,J,K)
	2792	CALL wrf_error_fatal(message)
	2793	ENDIF
	2794
	2795	ENDDO
	2796	ENDDO
	2797	ENDDO
	2798
	2799	END SUBROUTINE compute_3d_pressure
	2800
	2801	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2802	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2803	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2804
	2805	SUBROUTINE interp_press2press_lin(press_in,press_out, &
	2806	data_in, data_out,generic &
	2807	&, extrapolate,ignore_lowest,TFIELD &
	2808	&, IDS,IDE,JDS,JDE,KDS,KDE &
	2809	&, IMS,IME,JMS,JME,KMS,KME &
	2810	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time )
	2811
	2812	! Interpolates data from one set of pressure surfaces to
	2813	! another set of pressures
	2814
	2815	INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
	2816	INTEGER :: IMS,IME,JMS,JME,KMS,KME
	2817	INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
	2818	INTEGER :: internal_time
	2819
	2820	! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME)
	2821	REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
	2822	REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
	2823	! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME)
	2824	REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic)
	2825	REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME)
	2826	LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest, TFIELD
	2827	LOGICAL :: col_smooth
	2828
	2829	INTEGER :: i,j
	2830	INTEGER :: k,kk
	2831	REAL :: desired_press
	2832	REAL :: dvaldlnp,dlnp,tadiabat,tiso
	2833
	2834	REAL, PARAMETER :: ADIAFAC=9.81/1004.
	2835	REAL, PARAMETER :: TSTEXTRAPFAC=.0065
	2836
	2837
	2838
	2839	DO K=KMS,KME
	2840	DO J=JMS,JME
	2841	DO I=IMS,IME
	2842	DATA_OUT(I,J,K)=-99999.9
	2843	ENDDO
	2844	ENDDO
	2845	ENDDO
	2846
	2847	IF (ignore_lowest) then
	2848	LMIN=2
	2849	ELSE
	2850	LMIN=1
	2851	ENDIF
	2852
	2853	DO j = JTS, min(JTE,JDE-1)
	2854	test_i: DO i = ITS, min(ITE,IDE-1)
	2855
	2856	IF (internal_time_loop .gt. 1) THEN
	2857	IF (J .ne. JDS .and. J .ne. JDE-1 .and. &
	2858	I .ne. IDS .and. I .ne. IDE-1 ) THEN
	2859	!! not on boundary
	2860	CYCLE test_i
	2861	ENDIF
	2862	ENDIF
	2863
	2864
	2865	col_smooth=.false.
	2866
	2867	output_loop: DO k = KDS,KDE-1
	2868
	2869	desired_press = press_out(i,j,k)
	2870
	2871	if (K .gt. KDS) then
	2872	if (TFIELD .and. col_smooth .and. desired_press .le. press_in(i,j,LMIN) &
	2873	.and. press_out(i,j,k-1) .ge. press_in(i,j,LMIN)) then
	2874	MAX_SMOOTH=K
	2875	! write(message,*) 'I,J, MAX_SMOOTH: ', I,J, MAX_SMOOTH
	2876	! CALL wrf_debug(100,message)
	2877	endif
	2878	endif
	2879
	2880	! keep track of where the extrapolation begins
	2881
	2882	IF (desired_press .GT. press_in(i,j,LMIN)) THEN
	2883	IF (TFIELD .and. K .eq. 1 .and. (desired_press - press_in(i,j,LMIN)) .gt. 3000.) then
	2884	col_smooth=.TRUE. ! due to large extrapolation distance
	2885	ENDIF
	2886
	2887
	2888	IF ((desired_press - press_in(i,j,LMIN)).LT. 50.) THEN ! 0.5 mb
	2889	data_out(i,j,k) = data_in(i,j,LMIN)
	2890	ELSE
	2891	IF (extrapolate) THEN
	2892	! Extrapolate downward because desired P level is below
	2893	! the lowest level in our input data. Extrapolate using simple
	2894	! 1st derivative of value with respect to ln P for the bottom 2
	2895	! input layers.
	2896
	2897	! Add a check to make sure we are not using the gradient of
	2898	! a very thin layer
	2899
	2900	if (TFIELD) then
	2901	tiso=0.5*(data_in(i,j,1)+data_in(i,j,2))
	2902	endif
	2903
	2904
	2905	IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 500.) THEN ! likely isobaric data
	2906	dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1))
	2907	dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+1)) / dlnp
	2908	ELSE ! assume terrain following
	2909	dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+5))
	2910	dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+5)) / dlnp
	2911	ENDIF
	2912	data_out(i,j,k) = data_in(i,j,LMIN) + dvaldlnp * &
	2913	( log(desired_press)-log(press_in(i,j,LMIN)) )
	2914
	2915	if (TFIELD .and. data_out(i,j,k) .lt. tiso-0.2) then
	2916
	2917	! restrict slope to -1K/10 hPa
	2918	dvaldlnp=max(dvaldlnp, -1.0/ &
	2919	log( press_in(i,j,LMIN) / &
	2920	( press_in(i,j,LMIN)-1000.) ))
	2921
	2922	data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * &
	2923	( log(desired_press)-log(press_in(i,j,LMIN)) )
	2924
	2925	elseif (TFIELD .and. data_out(i,j,k) .gt. tiso+0.2) then
	2926
	2927	! restrict slope to +0.8K/10 hPa
	2928	dvaldlnp=min(dvaldlnp, 0.8/ &
	2929	log( press_in(i,j,LMIN) / &
	2930	( press_in(i,j,LMIN)-1000.) ))
	2931
	2932	data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * &
	2933	( log(desired_press)-log(press_in(i,j,LMIN)) )
	2934
	2935	endif
	2936
	2937	ELSE
	2938	data_out(i,j,k) = data_in(i,j,LMIN)
	2939	ENDIF
	2940	ENDIF
	2941	ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN
	2942	IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN
	2943	data_out(i,j,k) = data_in(i,j,generic)
	2944	ELSE
	2945	IF (extrapolate) THEN
	2946	! Extrapolate upward
	2947	IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN
	2948	dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1))
	2949	dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-1))/dlnp
	2950	ELSE
	2951	dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2))
	2952	dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-2))/dlnp
	2953	ENDIF
	2954	data_out(i,j,k) = data_in(i,j,generic) + &
	2955	dvaldlnp * (log(desired_press)-log(press_in(i,j,generic)))
	2956	ELSE
	2957	data_out(i,j,k) = data_in(i,j,generic)
	2958	ENDIF
	2959	ENDIF
	2960	ELSE
	2961	! We can trap between two levels and linearly interpolate
	2962
	2963	input_loop: DO kk = LMIN, generic-1
	2964	IF (desired_press .EQ. press_in(i,j,kk) )THEN
	2965	data_out(i,j,k) = data_in(i,j,kk)
	2966	EXIT input_loop
	2967	ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. &
	2968	(desired_press .GT. press_in(i,j,kk+1)) ) THEN
	2969
	2970	! do trapped in lnp
	2971
	2972	dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1))
	2973	dvaldlnp = (data_in(i,j,kk)-data_in(i,j,kk+1))/dlnp
	2974	data_out(i,j,k) = data_in(i,j,kk+1)+ &
	2975	dvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1)))
	2976
	2977	EXIT input_loop
	2978	ENDIF
	2979
	2980	ENDDO input_loop
	2981	ENDIF
	2982	ENDDO output_loop
	2983
	2984	if (col_smooth) then
	2985	do K=max(KDS,MAX_SMOOTH-4),MAX_SMOOTH+4
	2986	data_out(I,J,K)=0.5*(data_out(I,J,K)+data_out(I,J,K+1))
	2987	enddo
	2988	endif
	2989
	2990	ENDDO test_i
	2991	ENDDO
	2992	END SUBROUTINE interp_press2press_lin
	2993
	2994	SUBROUTINE wind_adjust(press_in,press_out, &
	2995	U_in, V_in,U_out,V_out &
	2996	&, generic,depth_replace &
	2997	&, IDS,IDE,JDS,JDE,KDS,KDE &
	2998	&, IMS,IME,JMS,JME,KMS,KME &
	2999	&, ITS,ITE,JTS,JTE,KTS,KTE )
	3000
	3001	INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
	3002	INTEGER :: IMS,IME,JMS,JME,KMS,KME
	3003	INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
	3004	INTEGER :: MAXLIN,MAXLOUT
	3005
	3006	REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
	3007	REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
	3008	REAL, INTENT(IN) :: U_in(IMS:IME,JMS:JME,generic)
	3009	REAL, INTENT(IN) :: V_in(IMS:IME,JMS:JME,generic)
	3010	REAL, INTENT(INOUT) :: U_out(IMS:IME,KMS:KME,JMS:JME)
	3011	REAL, INTENT(INOUT) :: V_out(IMS:IME,KMS:KME,JMS:JME)
	3012	REAL :: p1d_in(generic)
	3013	REAL :: p1d_out(KDS:KDE-1)
	3014
	3015
	3016	DO j = JTS, min(JTE,JDE-1)
	3017	DO i = ITS, min(ITE,IDE-1)
	3018
	3019	! IF (press_out(I,J,1) .lt. press_in(I,J,2)) then
	3020	IF( (press_in(I,J,2)-press_out(I,J,1)) .gt. 200.) then
	3021
	3022	U_out(I,1,J)=U_in(I,J,2)
	3023	V_out(I,1,J)=V_in(I,J,2)
	3024
	3025	INLOOP: DO L=2,generic
	3026	p1d_in(L)=-9999.
	3027	IF ( (press_in(I,J,2)-press_in(I,J,L)) .lt. depth_replace) THEN
	3028	p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L))
	3029	MAXLIN=L
	3030	ELSE
	3031	p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L))
	3032	EXIT INLOOP
	3033	ENDIF
	3034	END DO INLOOP
	3035
	3036	OUTLOOP: DO L=KDS,KDE-1
	3037	p1d_out(L)=-9999.
	3038	IF ( (press_out(I,J,1)-press_out(I,J,L)) .lt. depth_replace) THEN
	3039	p1d_out(L)=(press_out(I,J,1)-press_out(I,J,L))
	3040	MAXLOUT=L
	3041	ELSE
	3042	EXIT OUTLOOP
	3043	ENDIF
	3044	END DO OUTLOOP
	3045
	3046	DO L=1,MAXLOUT
	3047	ptarg=p1d_out(L)
	3048
	3049	FINDLOOP: DO LL=2,MAXLIN
	3050
	3051	if (p1d_in(LL) .lt. ptarg .and. p1d_in(LL+1) .gt. ptarg) then
	3052
	3053	dlnp=log(p1d_in(LL))-log(p1d_in(LL+1))
	3054	dudlnp=(U_in(I,J,LL)-U_in(I,J,LL+1))/dlnp
	3055	dvdlnp=(V_in(I,J,LL)-V_in(I,J,LL+1))/dlnp
	3056	U_out(I,L,J)=U_in(I,J,LL)+dudlnp*(log(ptarg)-log(p1d_in(LL)))
	3057	V_out(I,L,J)=V_in(I,J,LL)+dvdlnp*(log(ptarg)-log(p1d_in(LL)))
	3058
	3059	EXIT FINDLOOP
	3060	endif
	3061
	3062	END DO FINDLOOP
	3063	END DO ! MAXLOUT loop
	3064
	3065
	3066	ENDIF
	3067
	3068	ENDDO
	3069	ENDDO
	3070
	3071
	3072
	3073	END SUBROUTINE wind_adjust
	3074	!--------------------------------------------------------------------
	3075
	3076	SUBROUTINE interp_press2press_log(press_in,press_out, &
	3077	data_in, data_out, generic &
	3078	&, extrapolate,ignore_lowest &
	3079	&, IDS,IDE,JDS,JDE,KDS,KDE &
	3080	&, IMS,IME,JMS,JME,KMS,KME &
	3081	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time )
	3082
	3083	! Interpolates ln(data) from one set of pressure surfaces to
	3084	! another set of pressures
	3085
	3086	INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
	3087	INTEGER :: IMS,IME,JMS,JME,KMS,KME
	3088	INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
	3089	INTEGER :: internal_time
	3090
	3091	! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME)
	3092	REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
	3093	REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
	3094	! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME)
	3095	! REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic)
	3096	REAL :: data_in(IMS:IME,JMS:JME,generic)
	3097	REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME)
	3098	LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest
	3099
	3100	INTEGER :: i,j
	3101	INTEGER :: k,kk
	3102	REAL :: desired_press
	3103	REAL :: dlnvaldlnp,dlnp
	3104
	3105
	3106	DO K=1,generic
	3107	DO j = JTS, min(JTE,JDE-1)
	3108	DO i = ITS, min(ITE,IDE-1)
	3109	DATA_IN(I,J,K)=max(DATA_in(I,J,K),1.e-13)
	3110	ENDDO
	3111	ENDDO
	3112	ENDDO
	3113
	3114	DO K=KMS,KME
	3115	DO J=JMS,JME
	3116	DO I=IMS,IME
	3117	DATA_OUT(I,J,K)=-99999.9
	3118	ENDDO
	3119	ENDDO
	3120	ENDDO
	3121
	3122	IF (ignore_lowest) then
	3123	LMIN=2
	3124	ELSE
	3125	LMIN=1
	3126	ENDIF
	3127
	3128	DO j = JTS, min(JTE,JDE-1)
	3129	test_i: DO i = ITS, min(ITE,IDE-1)
	3130
	3131	IF (internal_time .gt. 1) THEN
	3132	IF (J .ne. JDS .and. J .ne. JDE-1 .and. &
	3133	I .ne. IDS .and. I .ne. IDE-1 ) THEN
	3134	!! not on boundary
	3135	CYCLE test_i
	3136	ENDIF
	3137	ENDIF
	3138
	3139
	3140	output_loop: DO k = KDS,KDE-1
	3141
	3142	desired_press = press_out(i,j,k)
	3143
	3144	IF (desired_press .GT. press_in(i,j,LMIN)) THEN
	3145
	3146	IF ((desired_press - press_in(i,j,LMIN)).LT. 10.) THEN ! 0.1 mb
	3147	data_out(i,j,k) = data_in(i,j,LMIN)
	3148	ELSE
	3149	IF (extrapolate) THEN
	3150	! Extrapolate downward because desired P level is below
	3151	! the lowest level in our input data. Extrapolate using simple
	3152	! 1st derivative of value with respect to ln P for the bottom 2
	3153	! input layers.
	3154
	3155	! Add a check to make sure we are not using the gradient of
	3156	! a very thin layer
	3157
	3158	IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 100.) THEN
	3159	dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1))
	3160	dlnvaldlnp = ( log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+1)) ) / dlnp
	3161
	3162	ELSE
	3163
	3164	dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+2))
	3165	dlnvaldlnp = (log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+2))) / dlnp
	3166
	3167	ENDIF
	3168
	3169	data_out(i,j,k) = exp(log(data_in(i,j,LMIN)) + dlnvaldlnp * &
	3170	( log(desired_press)-log(press_in(i,j,LMIN))))
	3171	ELSE
	3172	data_out(i,j,k) = data_in(i,j,LMIN)
	3173	ENDIF
	3174	ENDIF
	3175	ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN
	3176	IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN
	3177	data_out(i,j,k) = data_in(i,j,generic)
	3178	ELSE
	3179	IF (extrapolate) THEN
	3180	! Extrapolate upward
	3181	IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN
	3182	dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1))
	3183	dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-1)))/dlnp
	3184	ELSE
	3185	dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2))
	3186	dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-2)))/dlnp
	3187	ENDIF
	3188	data_out(i,j,k) = exp(log(data_in(i,j,generic)) + &
	3189	dlnvaldlnp * (log(desired_press)-log(press_in(i,j,generic))))
	3190	ELSE
	3191	data_out(i,j,k) = data_in(i,j,generic)
	3192	ENDIF
	3193	ENDIF
	3194	ELSE
	3195	! We can trap between two levels and linearly interpolate
	3196
	3197	input_loop: DO kk = LMIN, generic-1
	3198	IF (desired_press .EQ. press_in(i,j,kk) )THEN
	3199	data_out(i,j,k) = data_in(i,j,kk)
	3200	EXIT input_loop
	3201	ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. &
	3202	(desired_press .GT. press_in(i,j,kk+1)) ) THEN
	3203
	3204	! do trapped in lnp
	3205
	3206	dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1))
	3207	dlnvaldlnp = (log(data_in(i,j,kk))-log(data_in(i,j,kk+1)))/dlnp
	3208	data_out(i,j,k) = exp(log(data_in(i,j,kk+1))+ &
	3209	dlnvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1))))
	3210
	3211	EXIT input_loop
	3212
	3213	ENDIF
	3214
	3215	ENDDO input_loop
	3216	ENDIF
	3217	ENDDO output_loop
	3218	ENDDO test_i
	3219	ENDDO
	3220	END SUBROUTINE interp_press2press_log
	3221
	3222	!-------------------------------------------------------------------
	3223	SUBROUTINE rh_to_mxrat (rh, t, p, q , wrt_liquid , &
	3224	ids , ide , jds , jde , kds , kde , &
	3225	ims , ime , jms , jme , kms , kme , &
	3226	its , ite , jts , jte , kts , kte )
	3227
	3228	IMPLICIT NONE
	3229
	3230	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
	3231	ims , ime , jms , jme , kms , kme , &
	3232	its , ite , jts , jte , kts , kte
	3233
	3234	LOGICAL , INTENT(IN) :: wrt_liquid
	3235
	3236	! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
	3237	! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
	3238	REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(IN) :: p , t
	3239	REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(INOUT) :: rh
	3240	REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(OUT) :: q
	3241
	3242	! Local vars
	3243
	3244	INTEGER :: i , j , k
	3245
	3246	REAL :: ew , q1 , t1
	3247
	3248	REAL, PARAMETER :: T_REF = 0.0
	3249	REAL, PARAMETER :: MW_AIR = 28.966
	3250	REAL, PARAMETER :: MW_VAP = 18.0152
	3251
	3252	REAL, PARAMETER :: A0 = 6.107799961
	3253	REAL, PARAMETER :: A1 = 4.436518521e-01
	3254	REAL, PARAMETER :: A2 = 1.428945805e-02
	3255	REAL, PARAMETER :: A3 = 2.650648471e-04
	3256	REAL, PARAMETER :: A4 = 3.031240396e-06
	3257	REAL, PARAMETER :: A5 = 2.034080948e-08
	3258	REAL, PARAMETER :: A6 = 6.136820929e-11
	3259
	3260	REAL, PARAMETER :: ES0 = 6.1121
	3261
	3262	REAL, PARAMETER :: C1 = 9.09718
	3263	REAL, PARAMETER :: C2 = 3.56654
	3264	REAL, PARAMETER :: C3 = 0.876793
	3265	REAL, PARAMETER :: EIS = 6.1071
	3266	REAL :: RHS
	3267	REAL, PARAMETER :: TF = 273.16
	3268	REAL :: TK
	3269
	3270	REAL :: ES
	3271	REAL :: QS
	3272	REAL, PARAMETER :: EPS = 0.622
	3273	REAL, PARAMETER :: SVP1 = 0.6112
	3274	REAL, PARAMETER :: SVP2 = 17.67
	3275	REAL, PARAMETER :: SVP3 = 29.65
	3276	REAL, PARAMETER :: SVPT0 = 273.15
	3277
	3278	! This subroutine computes mixing ratio (q, kg/kg) from basic variables
	3279	! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%).
	3280	! The reference temperature (t_ref, C) is used to describe the temperature
	3281	! at which the liquid and ice phase change occurs.
	3282
	3283	DO k = kts , kte
	3284	DO j = jts , MIN ( jde-1 , jte )
	3285	DO i = its , MIN (ide-1 , ite )
	3286	rh(i,j,k) = MIN ( MAX ( rh(i,j,k) , 1. ) , 100. )
	3287	END DO
	3288	END DO
	3289	END DO
	3290
	3291	IF ( wrt_liquid ) THEN
	3292	DO k = kts , kte
	3293	DO j = jts , MIN ( jde-1 , jte )
	3294	DO i = its , MIN (ide-1 , ite )
	3295	es=svp110.EXP(svp2*(t(i,j,k)-svpt0)/(t(i,j,k)-svp3))
	3296	qs=eps*es/(p(i,j,k)/100.-es)
	3297	q(i,j,k)=MAX(.01rh(i,j,k)qs,0.0)
	3298	END DO
	3299	END DO
	3300	END DO
	3301
	3302	ELSE
	3303	DO k = kts , kte
	3304	DO j = jts , MIN ( jde-1 , jte )
	3305	DO i = its , MIN (ide-1 , ite )
	3306
	3307	t1 = t(i,j,k) - 273.16
	3308
	3309	! Obviously dry.
	3310
	3311	IF ( t1 .lt. -200. ) THEN
	3312	q(i,j,k) = 0
	3313
	3314	ELSE
	3315
	3316	! First compute the ambient vapor pressure of water
	3317
	3318	IF ( ( t1 .GE. t_ref ) .AND. ( t1 .GE. -47.) ) THEN ! liq phase ESLO
	3319	ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6)))))
	3320
	3321	ELSE IF ( ( t1 .GE. t_ref ) .AND. ( t1 .LT. -47. ) ) then !liq phas poor ES
	3322	ew = es0 * exp(17.67 * t1 / ( t1 + 243.5))
	3323
	3324	ELSE
	3325	tk = t(i,j,k)
	3326	rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + &
	3327	c3 * (1. - tk / tf) + alog10(eis)
	3328	ew = 10. ** rhs
	3329
	3330	END IF
	3331
	3332	! Now sat vap pres obtained compute local vapor pressure
	3333
	3334	ew = MAX ( ew , 0. ) * rh(i,j,k) * 0.01
	3335
	3336	! Now compute the specific humidity using the partial vapor
	3337	! pressures of water vapor (ew) and dry air (p-ew). The
	3338	! constants assume that the pressure is in hPa, so we divide
	3339	! the pressures by 100.
	3340
	3341	q1 = mw_vap * ew
	3342	q1 = q1 / (q1 + mw_air * (p(i,j,k)/100. - ew))
	3343
	3344	q(i,j,k) = q1 / (1. - q1 )
	3345
	3346	END IF
	3347
	3348	END DO
	3349	END DO
	3350	END DO
	3351
	3352	END IF
	3353
	3354	END SUBROUTINE rh_to_mxrat
	3355
	3356	!--=------------------------------------------------------------------
	3357
	3358	SUBROUTINE boundary_smooth(h, landmask, grid, nsmth , nrow &
	3359	&, IDS,IDE,JDS,JDE,KDS,KDE &
	3360	&, IMS,IME,JMS,JME,KMS,KME &
	3361	&, ITS,ITE,JTS,JTE,KTS,KTE )
	3362
	3363	implicit none
	3364
	3365	TYPE (domain) :: grid
	3366
	3367	integer :: IDS,IDE,JDS,JDE,KDS,KDE
	3368	integer :: IMS,IME,JMS,JME,KMS,KME
	3369	integer :: ITS,ITE,JTS,JTE,KTS,KTE
	3370	integer:: ihw(JDS:JDE-1),ihe(JDS:JDE-1),nsmth,nrow
	3371	real:: h(IMS:IME,JMS:JME),landmask(IMS:IME,JMS:JME)
	3372	real :: h_old(IMS:IME,JMS:JME)
	3373	real :: hbms(IDS:IDE-1,JDS:JDE-1)
	3374	real :: hse(IDS:IDE-1,JDS:JDE-1)
	3375	real :: hne(IDS:IDE-1,JDS:JDE-1)
	3376	integer :: IPS,IPE,JPS,JPE,KPS,KPE
	3377	integer :: ihl, ihh, m2l, ibas,jmelin
	3378	integer :: I,J,KS,IOFFSET,JSTART,JEND
	3379	character (len=255) :: message
	3380
	3381	ips=its
	3382	ipe=ite
	3383	jps=jts
	3384	jpe=jte
	3385	kps=kts
	3386	kpe=kte
	3387
	3388	do j= JTS,min(JTE,JDE-1)
	3389	ihw(J)=-mod(J,2)
	3390	ihe(j)=ihw(J)+1
	3391	end do
	3392
	3393	do J=JTS,min(JTE,JDE-1)
	3394	do I=ITS,min(ITE,IDE-1)
	3395	hbms(I,J)=landmask(I,J)
	3396	enddo
	3397	enddo
	3398
	3399	jmelin=(JDE-1)-nrow+1
	3400	ibas=nrow/2
	3401	m2l=mod(nrow,2)
	3402
	3403	do j=jts,min(jte,jde-1)
	3404	ihl=ibas+mod(j,2)+m2l*mod(J+1,2)
	3405	ihh=(IDE-1)-ibas-m2l*mod(J+1,2)
	3406	do i=its,min(ite,ide-1)
	3407	if (I .ge. ihl .and. I .le. ihh .and. J .ge. nrow .and. J .le. jmelin) then
	3408	hbms(I,J)=0.
	3409	endif
	3410	end do
	3411	end do
	3412
	3413	634 format(30(f2.0,1x))
	3414
	3415	do KS=1,nsmth
	3416
	3417	grid%ht_gc=h
	3418	#ifdef DM_PARALLEL
	3419	# include "HALO_NMM_MG.inc"
	3420	#endif
	3421	h=grid%ht_gc
	3422	h_old=grid%ht_gc
	3423
	3424	do J=JTS,min(JTE,JDE-1)
	3425	do I=ITS, min(ITE,IDE-1)
	3426	if (I .ge. (IDS+mod(J,2)) .and. J .gt. JDS .and. J .lt. JDE-1 .and. I .lt. IDE-1) then
	3427	h(i,j)= ( h_old(i+ihe(j),j+1) + h_old(i+ihw(j),j-1) + h_old(i+ihe(j),j-1) + h_old(i+ihw(j),j+1) - &
	3428	4. h_old(i,j) )hbms(i,j)*0.125+h_old(i,j)
	3429	endif
	3430
	3431	enddo
	3432	enddo
	3433
	3434	! special treatment for four corners
	3435
	3436	if (hbms(1,1) .eq. 1 .and. ITS .le. 1 .and. JTS .le. 1) then
	3437	h(1,1)=0.75h(1,1)+0.125h(1+ihe(1),2)+ &
	3438	0.0625*(h(2,1)+h(1,3))
	3439	endif
	3440
	3441	if (hbms(IDE-1,1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTS .le. 1) then
	3442	h(IDE-1,1)=0.75h(IDE-1,1)+0.125h(IDE-1+ihw(1),2)+ &
	3443	0.0625*(h(IDE-1-1,1)+h(IDE-1,3))
	3444	endif
	3445
	3446	if (hbms(1,JDE-1) .eq. 1 .and. ITS .le. 1 .and. JTE .ge. JDE-2) then
	3447	h(1,JDE-1)=0.75h(1,JDE-1)+0.125h(1+ihe(JDE-1),JDE-1-1)+ &
	3448	0.0625*(h(2,JDE-1)+h(1,JDE-1-2))
	3449	endif
	3450
	3451	if (hbms(IDE-1,JDE-1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTE .ge. JDE-2) then
	3452	h(IDE-1,JDE-1)=0.75h(IDE-1,JDE-1)+0.125h(IDE-1+ihw(JDE-1),JDE-1-1)+ &
	3453	0.0625*(h(IDE-1-1,JDE-1)+h(IDE-1,JDE-1-2))
	3454	endif
	3455
	3456	do J=JMS,JME
	3457	do I=IMS,IME
	3458	grid%ht_gc(I,J)=h(I,J)
	3459	enddo
	3460	enddo
	3461	#ifdef DM_PARALLEL
	3462	# include "HALO_NMM_MG.inc"
	3463	#endif
	3464	do J=JMS,JME
	3465	do I=IMS,IME
	3466	h(I,J)=grid%ht_gc(I,J)
	3467	enddo
	3468	enddo
	3469
	3470
	3471	! S bound
	3472	if (JTS .eq. JDS) then
	3473	J=JTS
	3474
	3475	do I=ITS,ITE
	3476	if (I .ge. IDS+1 .and. I .le. IDE-2) then
	3477	if (hbms(I,J) .eq. 1) then
	3478	h(I,J)=0.75h(I,J)+0.125(h(I+ihw(J),J+1)+h(I+ihe(J),J+1))
	3479	endif
	3480	endif
	3481	enddo
	3482
	3483	endif
	3484
	3485	! N bound
	3486	if (JTE .eq. JDE) then
	3487	J=JDE-1
	3488	write(message,*) 'DOING N BOUND SMOOTHING for J= ', J
	3489	CALL wrf_message(message)
	3490	do I=ITS,min(ITE,IDE-1)
	3491	if (hbms(I,J) .eq. 1 .and. I .ge. IDS+1 .and. I .le. IDE-2) then
	3492	h(I,J)=0.75h(I,J)+0.125(h(I+ihw(J),J-1)+h(I+ihe(J),J-1))
	3493	endif
	3494	enddo
	3495	endif
	3496
	3497	! W bound
	3498	if (ITS .eq. IDS) then
	3499	I=ITS
	3500	do J=JTS,min(JTE,JDE-1)
	3501	if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then
	3502	h(I,J)=0.75h(I,J)+0.125(h(I+ihe(J),J+1)+h(I+ihe(J),J-1))
	3503	endif
	3504	enddo
	3505	endif
	3506
	3507	! E bound
	3508	if (ITE .eq. IDE) then
	3509	write(message,*) 'DOING E BOUND SMOOTHING for I= ', min(ITE,IDE-1)
	3510	CALL wrf_message(message)
	3511	I=min(ITE,IDE-1)
	3512	do J=JTS,min(JTE,JDE-1)
	3513	if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then
	3514	h(I,J)=0.75h(I,J)+0.125(h(I+ihw(J),J+1)+h(I+ihw(J),J-1))
	3515	endif
	3516	enddo
	3517	endif
	3518
	3519	enddo ! end ks loop
	3520
	3521	do J=JMS,JME
	3522	do I=IMS,IME
	3523	grid%ht_gc(I,J)=h(I,J)
	3524	enddo
	3525	enddo
	3526	#ifdef DM_PARALLEL
	3527	#include "HALO_NMM_MG.inc"
	3528	#endif
	3529	do J=JMS,JME
	3530	do I=IMS,IME
	3531	h(I,J)=grid%ht_gc(I,J)
	3532	enddo
	3533	enddo
	3534
	3535	! extra smoothing along inner boundary
	3536
	3537	if (JTS .eq. JDS) then
	3538	if (ITE .eq. IDE) then
	3539	IOFFSET=1
	3540	else
	3541	IOFFSET=0
	3542	endif
	3543	! Southern Boundary
	3544	do i=its,min(ITE,IDE-1)-IOFFSET
	3545	h(i,2)=0.25*(h(i,1)+h(i+1,1)+ &
	3546	h(i,3)+h(i+1,3))
	3547	enddo
	3548	endif
	3549
	3550
	3551	if (JTE .eq. JDE) then
	3552	if (ITE .eq. IDE) then
	3553	IOFFSET=1
	3554	else
	3555	IOFFSET=0
	3556	endif
	3557	do i=its,min(ITE,IDE-1)-IOFFSET
	3558	h(i,(JDE-1)-1)=0.25*(h(i,(JDE-1)-2)+h(i+1,(JDE-1)-2)+ &
	3559	h(i,JDE-1)+h(i+1,JDE-1))
	3560	enddo
	3561	endif
	3562
	3563	if (JTS .eq. 1) then
	3564	JSTART=4
	3565	else
	3566	JSTART=JTS+mod(JTS,2) ! needs to be even
	3567	endif
	3568
	3569	if (JTE .eq. JDE) then
	3570	JEND=(JDE-1)-3
	3571	else
	3572	JEND=JTE
	3573	endif
	3574
	3575	if (ITS .eq. IDS) then
	3576
	3577	! Western Boundary
	3578	do j=JSTART,JEND,2
	3579	h(1,j)=0.25*(h(1,j-1)+h(2,j-1)+ &
	3580	h(1,j+1)+h(2,j+1))
	3581
	3582	enddo
	3583	endif
	3584
	3585
	3586	if (ITE .eq. IDE) then
	3587	! Eastern Boundary
	3588	do j=JSTART,JEND,2
	3589	h((IDE-1)-1,j)=0.25*(h((IDE-1)-1,j-1)+h((IDE-1),j-1)+ &
	3590	h((IDE-1)-1,j+1)+h((IDE-1),j+1))
	3591	enddo
	3592	endif
	3593
	3594
	3595	END SUBROUTINE boundary_smooth
	3596
	3597	!--------------------------------------------------------------------
	3598
	3599	SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , &
	3600	ids , ide , jds , jde , kds , kde , &
	3601	ims , ime , jms , jme , kms , kme , &
	3602	its , ite , jts , jte , kts , kte )
	3603
	3604	! Linrarly in time interpolate data to a current valid time. The data is
	3605	! assumed to come in "monthly", valid at the 15th of every month.
	3606
	3607	IMPLICIT NONE
	3608
	3609	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
	3610	ims , ime , jms , jme , kms , kme , &
	3611	its , ite , jts , jte , kts , kte
	3612
	3613	CHARACTER (LEN=24) , INTENT(IN) :: date_str
	3614	REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in
	3615	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out
	3616
	3617	! Local vars
	3618
	3619	INTEGER :: i , j , l
	3620	INTEGER , DIMENSION(0:13) :: middle
	3621	INTEGER :: target_julyr , target_julday , target_date
	3622	INTEGER :: julyr , julday , int_month, next_month
	3623	REAL :: gmt
	3624	CHARACTER (LEN=4) :: yr
	3625	CHARACTER (LEN=2) :: mon , day15
	3626
	3627
	3628	WRITE(day15,FMT='(I2.2)') 15
	3629	DO l = 1 , 12
	3630	WRITE(mon,FMT='(I2.2)') l
	3631	CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt )
	3632	middle(l) = julyr*1000 + julday
	3633	END DO
	3634
	3635	l = 0
	3636	middle(l) = middle( 1) - 31
	3637
	3638	l = 13
	3639	middle(l) = middle(12) + 31
	3640
	3641	CALL get_julgmt ( date_str , target_julyr , target_julday , gmt )
	3642	target_date = target_julyr * 1000 + target_julday
	3643	find_month : DO l = 0 , 12
	3644	IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN
	3645	DO j = jts , MIN ( jde-1 , jte )
	3646	DO i = its , MIN (ide-1 , ite )
	3647	int_month = MOD ( l , 12 )
	3648	IF ( int_month .EQ. 0 ) int_month = 12
	3649
	3650	IF (int_month == 12) THEN
	3651	next_month=1
	3652	ELSE
	3653	next_month=int_month+1
	3654	ENDIF
	3655
	3656	field_out(i,j) = ( field_in(i,j,next_month) * ( target_date - middle(l) ) + &
	3657	field_in(i,j,int_month ) * ( middle(l+1) - target_date ) ) / &
	3658	( middle(l+1) - middle(l) )
	3659	END DO
	3660	END DO
	3661	EXIT find_month
	3662	END IF
	3663	END DO find_month
	3664	END SUBROUTINE monthly_interp_to_date
	3665
	3666	!---------------------------------------------------------------------
	3667	SUBROUTINE monthly_min_max ( field_in , field_min , field_max , &
	3668	ids , ide , jds , jde , kds , kde , &
	3669	ims , ime , jms , jme , kms , kme , &
	3670	its , ite , jts , jte , kts , kte )
	3671
	3672	! Plow through each month, find the max, min values for each i,j.
	3673
	3674	IMPLICIT NONE
	3675
	3676	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
	3677	ims , ime , jms , jme , kms , kme , &
	3678	its , ite , jts , jte , kts , kte
	3679
	3680	REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in
	3681	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max
	3682
	3683	! Local vars
	3684
	3685	INTEGER :: i , j , l
	3686	REAL :: minner , maxxer
	3687
	3688	DO j = jts , MIN(jde-1,jte)
	3689	DO i = its , MIN(ide-1,ite)
	3690	minner = field_in(i,j,1)
	3691	maxxer = field_in(i,j,1)
	3692	DO l = 2 , 12
	3693	IF ( field_in(i,j,l) .LT. minner ) THEN
	3694	minner = field_in(i,j,l)
	3695	END IF
	3696	IF ( field_in(i,j,l) .GT. maxxer ) THEN
	3697	maxxer = field_in(i,j,l)
	3698	END IF
	3699	END DO
	3700	field_min(i,j) = minner
	3701	field_max(i,j) = maxxer
	3702	END DO
	3703	END DO
	3704
	3705	END SUBROUTINE monthly_min_max
	3706
	3707	!-----------------------------------------------------------------------
	3708
	3709	SUBROUTINE reverse_vert_coord ( field, start_z, end_z &
	3710	&, IDS,IDE,JDS,JDE,KDS,KDE &
	3711	&, IMS,IME,JMS,JME,KMS,KME &
	3712	&, ITS,ITE,JTS,JTE,KTS,KTE )
	3713
	3714	IMPLICIT NONE
	3715
	3716	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
	3717	ims , ime , jms , jme , kms , kme , &
	3718	its , ite , jts , jte , kts , kte, &
	3719	start_z, end_z
	3720
	3721	REAL, INTENT(INOUT) :: field(IMS:IME,JMS:JME,end_z)
	3722	! local
	3723
	3724	INTEGER :: I,J,L
	3725	REAL, ALLOCATABLE :: dum3d(:,:,:)
	3726
	3727	allocate(dum3d(IMS:IME,JMS:JME,end_z))
	3728
	3729	DO L=start_z,end_z
	3730	DO J=jts,min(jte,jde-1)
	3731	DO I=its,min(ite,ide-1)
	3732	dum3d(I,J,L)=field(I,J,end_z-L+start_z)
	3733	END DO
	3734	END DO
	3735	END DO
	3736
	3737	DO L=start_z,end_z
	3738	DO J=jts,min(jte,jde-1)
	3739	DO I=its,min(ite,ide-1)
	3740	field(I,J,L)=dum3d(I,J,L)
	3741	END DO
	3742	END DO
	3743	END DO
	3744
	3745	DEALLOCATE(dum3d)
	3746
	3747	END SUBROUTINE reverse_vert_coord
	3748
	3749
	3750	!--------------------------------------------------------------------
	3751
	3752	SUBROUTINE compute_nmm_levels(ninterface, ptop, eta_levels)
	3753
	3754	USE module_model_constants
	3755
	3756	IMPLICIT NONE
	3757
	3758	character(len=132):: message
	3759	integer :: ninterface,Lthick,L
	3760	real, parameter:: gamma=.0065
	3761	real, parameter:: t_stand=288.
	3762	real, parameter:: p_stand=101325.
	3763
	3764	real :: maxdz_compute, ptop
	3765	real :: plower,pupper,tlay, sum
	3766
	3767	real :: eta_levels(ninterface)
	3768	real, allocatable:: Z(:)
	3769	real, allocatable:: deta_levels_spline(:)
	3770
	3771	logical:: print_pbl_warn
	3772
	3773	!----------------------------------------------------
	3774
	3775	allocate(Z(ninterface))
	3776	allocate(deta_levels_spline(ninterface-1))
	3777
	3778	CALL compute_eta_spline(ninterface-1,deta_levels_spline,ptop)
	3779
	3780	sum=0.
	3781	DO L=1,ninterface-1
	3782	sum=sum+deta_levels_spline(L)
	3783	ENDDO
	3784
	3785	eta_levels(1)=1.00
	3786
	3787	DO L=2,ninterface
	3788	eta_levels(L)=eta_levels(L-1)-deta_levels_spline(L-1)
	3789	ENDDO
	3790
	3791	eta_levels(ninterface)=0.00
	3792
	3793	DO L=2,ninterface-1
	3794	eta_levels(L)=0.5(eta_levels(L))+0.25(eta_levels(L-1)+eta_levels(L+1))
	3795	ENDDO
	3796
	3797	Z(1)=0.
	3798	maxdz_compute=0.
	3799	print_pbl_warn=.false.
	3800
	3801	DO L=2,ninterface
	3802	tlay=max( t_stand-gamma*Z(L-1), 216.5)
	3803	plower=ptop+(p_stand-ptop)*eta_levels(L-1)
	3804	pupper=ptop+(p_stand-ptop)*eta_levels(L)
	3805	Z(L)=Z(L-1)+(tlayr_d/g)(log(plower)-log(pupper))
	3806
	3807	if (plower .gt. 85000. .and. pupper .lt. 85000. .and. L .lt. 10) then
	3808	print_pbl_warn=.true.
	3809	endif
	3810
	3811	write(message,*) 'L, eta(l), pupper, Z(L): ', L, eta_levels(L),pupper,Z(L)
	3812	CALL wrf_debug(100,message)
	3813
	3814	if (Z(L)-Z(L-1) .gt. maxdz_compute) then
	3815	Lthick=L
	3816	endif
	3817
	3818	maxdz_compute=max(maxdz_compute,Z(L)-Z(L-1))
	3819	ENDDO
	3820
	3821	if (print_pbl_warn) then
	3822	write(message,*) 'WARNING - PBL MAY BE POORLY RESOLVED WITH NUMBER OF VERTICAL LEVELS'
	3823	CALL wrf_message(message)
	3824	write(message,*) ' - CONSIDER INCREASING THE VERTICAL RESOLUTION'
	3825	CALL wrf_message(message)
	3826	endif
	3827
	3828	write(message,*) 'thickest layer was: ', maxdz_compute , 'meters thick at level: ', Lthick
	3829	CALL wrf_message(message)
	3830
	3831	END SUBROUTINE compute_nmm_levels
	3832
	3833	!---------------------------
	3834
	3835	SUBROUTINE compute_eta_spline(LM, dsg, ptop)
	3836
	3837	IMPLICIT NONE
	3838
	3839	real:: dsg(LM), ptop, sum, rsum
	3840	real, allocatable:: xold(:),dold(:)
	3841	real, allocatable:: xnew(:),sgm(:)
	3842	real, allocatable:: pps(:),qqs(:),y2s(:)
	3843	integer nlev,LM,L,KOLD
	3844
	3845	IF (LM .ge. 46) THEN
	3846	KOLD=9
	3847	allocate(xold(KOLD))
	3848	allocate(dold(KOLD))
	3849
	3850	xold(1)=.00
	3851	dold(1)=.006
	3852	xold(2)=.13
	3853	dold(2)=.009
	3854	xold(3)=.19
	3855	dold(3)=.012
	3856	xold(4)=.30
	3857	dold(4)=.036
	3858	xold(5)=.42
	3859	dold(5)=.041
	3860	xold(6)=.56
	3861	dold(6)=.040
	3862	xold(7)=.69
	3863	dold(7)=.018
	3864
	3865	if (ptop .ge. 2000.) then
	3866	xold(8)=.90
	3867	dold(8)=.012
	3868	xold(9)=1.0
	3869	dold(9)=.006
	3870	else
	3871	xold(8)=.90
	3872	dold(8)=.008
	3873	xold(9)=1.0
	3874	dold(9)=.003
	3875	endif
	3876
	3877	ELSE
	3878
	3879	KOLD=8
	3880	allocate(xold(KOLD))
	3881	allocate(dold(KOLD))
	3882
	3883	xold(1)=.00
	3884	dold(1)=.006
	3885	xold(2)=.18
	3886	dold(2)=.015
	3887	xold(3)=.32
	3888	dold(3)=.035
	3889	xold(4)=.50
	3890	dold(4)=.040
	3891	xold(5)=.68
	3892	dold(5)=.030
	3893	xold(6)=.75
	3894	dold(6)=.017
	3895	xold(7)=.85
	3896	dold(7)=.012
	3897
	3898	if (ptop .ge. 2000.) then
	3899	xold(8)=1.0
	3900	dold(8)=.013
	3901	else
	3902	xold(8)=1.0
	3903	dold(8)=.008
	3904	endif
	3905
	3906	ENDIF
	3907
	3908	allocate(xnew(lm))
	3909	allocate(sgm(lm))
	3910	allocate(pps(lm))
	3911	allocate(qqs(lm))
	3912	allocate(y2s(lm))
	3913
	3914	DO L=1,LM
	3915	xnew(l)=float(l-1)/float(lm-1)
	3916	ENDDO
	3917
	3918	y2s=0.
	3919
	3920	CALL spline(kold,kold,xold,dold,y2s,lm,xnew,dsg,pps,qqs)
	3921
	3922	sum=0.
	3923	DO l=1,lm
	3924	sum=sum+dsg(l)
	3925	ENDDO
	3926
	3927	rsum=1./sum
	3928	sgm(1)=0.
	3929
	3930	DO L=1,lm-1
	3931	dsg(l)=dsg(l)*rsum
	3932	sgm(l+1)=sgm(l)+dsg(l)
	3933	ENDDO
	3934	sgm(lm+1)=1.
	3935	dsg(lm)=sgm(lm+1)-sgm(lm)
	3936
	3937	END SUBROUTINE compute_eta_spline
	3938
	3939	! -------------------------------------------------------------------
	3940
	3941	subroutine spline(JTBX,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
	3942	! ********************************************************************
	3943	! * *
	3944	! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
	3945	! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
	3946	! * *
	3947	! * PROGRAMER Z. JANJIC *
	3948	! * *
	3949	! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
	3950	! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
	3951	! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
	3952	! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
	3953	! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
	3954	! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
	3955	! * SPECIFIED. *
	3956	! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
	3957	! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
	3958	! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
	3959	! * AND LE XOLD(NOLD). *
	3960	! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
	3961	! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
	3962	! * *
	3963	! ********************************************************************
	3964	!
	3965	! LOG:
	3966	!
	3967	! PYLE - June 2007 - eliminated use of GO TO statements.
	3968	!
	3969	!-----------------------------------------------------------------------
	3970	IMPLICIT NONE
	3971	!-----------------------------------------------------------------------
	3972	INTEGER,INTENT(IN) :: JTBX,NNEW,NOLD
	3973	REAL,DIMENSION(JTBX),INTENT(IN) :: XNEW,XOLD,YOLD
	3974	REAL,DIMENSION(JTBX),INTENT(INOUT) :: P,Q,Y2
	3975	REAL,DIMENSION(JTBX),INTENT(OUT) :: YNEW
	3976	!
	3977	INTEGER :: K,K1,K2,KOLD,NOLDM1
	3978	REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR &
	3979	& ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
	3980	!-----------------------------------------------------------------------
	3981	NOLDM1=NOLD-1
	3982
	3983	DXL=XOLD(2)-XOLD(1)
	3984	DXR=XOLD(3)-XOLD(2)
	3985	DYDXL=(YOLD(2)-YOLD(1))/DXL
	3986	DYDXR=(YOLD(3)-YOLD(2))/DXR
	3987	RTDXC=0.5/(DXL+DXR)
	3988
	3989	P(1)= RTDXC(6.(DYDXR-DYDXL)-DXL*Y2(1))
	3990	Q(1)=-RTDXC*DXR
	3991
	3992	first_loop: DO K=3,NOLD-1
	3993	IF(NOLD==3) exit first_loop ! should be impossible to execute
	3994	DXL=DXR
	3995	DYDXL=DYDXR
	3996	DXR=XOLD(K+1)-XOLD(K)
	3997	DYDXR=(YOLD(K+1)-YOLD(K))/DXR
	3998	DXC=DXL+DXR
	3999	DEN=1./(DXL*Q(K-2)+DXC+DXC)
	4000	P(K-1)= DEN(6.(DYDXR-DYDXL)-DXL*P(K-2))
	4001	Q(K-1)=-DEN*DXR
	4002	END DO first_loop
	4003
	4004	DO K=NOLDM1,2,-1
	4005	Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)
	4006	END DO
	4007
	4008	!-----------------------------------------------------------------------
	4009	second_loop: DO K1=1,NNEW+1
	4010	XK=XNEW(K1)
	4011	third_loop: DO K2=2,NOLD
	4012	IF(XOLD(K2)>XK)THEN
	4013	KOLD=K2-1
	4014	exit third_loop
	4015	ENDIF
	4016	END DO third_loop
	4017
	4018	IF (XOLD(K2) .le. XK) THEN
	4019	YNEW(K1)=YOLD(NOLD)
	4020	CYCLE second_loop
	4021	ENDIF
	4022
	4023	IF (K1 .eq. 1 .or. K .ne. KOLD) THEN
	4024	K=KOLD
	4025	Y2K=Y2(K)
	4026	Y2KP1=Y2(K+1)
	4027	DX=XOLD(K+1)-XOLD(K)
	4028	RDX=1./DX
	4029	AK=.1666667RDX(Y2KP1-Y2K)
	4030	BK=0.5*Y2K
	4031	CK=RDX(YOLD(K+1)-YOLD(K))-.1666667DX*(Y2KP1+Y2K+Y2K)
	4032	ENDIF
	4033
	4034	X=XK-XOLD(K)
	4035	XSQ=X*X
	4036	YNEW(K1)=AKXSQX+BKXSQ+CKX+YOLD(K)
	4037	END DO second_loop
	4038
	4039	END SUBROUTINE SPLINE
	4040
	4041	!--------------------------------------------------------------------
	4042	SUBROUTINE NMM_SH2O(IMS,IME,JMS,JME,ISTART,IM,JSTART,JM,&
	4043	NSOIL,ISLTPK, &
	4044	SM,SICE,STC,SMC,SH2O)
	4045
	4046	!! INTEGER, PARAMETER:: NSOTYP=9
	4047	! INTEGER, PARAMETER:: NSOTYP=16
	4048	INTEGER, PARAMETER:: NSOTYP=19 !!!!!!!!MAYBE???
	4049
	4050	REAL :: PSIS(NSOTYP),BETA(NSOTYP),SMCMAX(NSOTYP)
	4051	REAL :: STC(IMS:IME,NSOIL,JMS:JME), &
	4052	SMC(IMS:IME,NSOIL,JMS:JME)
	4053	REAL :: SH2O(IMS:IME,NSOIL,JMS:JME),SICE(IMS:IME,JMS:JME),&
	4054	SM(IMS:IME,JMS:JME)
	4055	REAL :: HLICE,GRAV,T0,BLIM
	4056	INTEGER :: ISLTPK(IMS:IME,JMS:JME)
	4057	CHARACTER(LEN=255) :: message
	4058
	4059	! Constants used in cold start SH2O initialization
	4060	DATA HLICE/3.335E5/,GRAV/9.81/,T0/273.15/
	4061	DATA BLIM/5.5/
	4062	! DATA PSIS /0.04,0.62,0.47,0.14,0.10,0.26,0.14,0.36,0.04/
	4063	! DATA BETA /4.26,8.72,11.55,4.74,10.73,8.17,6.77,5.25,4.26/
	4064	! DATA SMCMAX /0.421,0.464,0.468,0.434,0.406, &
	4065	! 0.465,0.404,0.439,0.421/
	4066
	4067
	4068	!!! NOT SURE...PSIS=SATPSI, BETA=BB??
	4069
	4070	DATA PSIS /0.069, 0.036, 0.141, 0.759, 0.759, 0.355, &
	4071	0.135, 0.617, 0.263, 0.098, 0.324, 0.468, &
	4072	0.355, 0.000, 0.069, 0.036, 0.468, 0.069, 0.069 /
	4073
	4074	DATA BETA/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, &
	4075	6.66, 8.72, 8.17, 10.73, 10.39, 11.55, &
	4076	5.25, 0.00, 2.79, 4.26, 11.55, 2.79, 2.79 /
	4077
	4078	DATA SMCMAX/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
	4079	0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
	4080	0.439, 1.000, 0.200, 0.421, 0.468, 0.200, 0.339/
	4081
	4082	DO K=1,NSOIL
	4083	DO J=JSTART,JM
	4084	DO I=ISTART,IM
	4085
	4086	!tst
	4087	IF (SMC(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
	4088	if (K .eq. 1) then
	4089	write(message,*) 'I,J,reducing SMC from ' ,I,J,SMC(I,K,J), 'to ', SMCMAX(ISLTPK(I,J))
	4090	CALL wrf_debug(100,message)
	4091	endif
	4092	SMC(I,K,J)=SMCMAX(ISLTPK(I,J))
	4093	ENDIF
	4094	!tst
	4095
	4096	IF ( (SM(I,J) .lt. 0.5) .and. (SICE(I,J) .lt. 0.5) ) THEN
	4097
	4098	IF (ISLTPK(I,J) .gt. 19) THEN
	4099	WRITE(message,*) 'FORCING ISLTPK at : ', I,J
	4100	CALL wrf_message(message)
	4101	ISLTPK(I,J)=9
	4102	ELSEIF (ISLTPK(I,J) .le. 0) then
	4103	WRITE(message,*) 'FORCING ISLTPK at : ', I,J
	4104	CALL wrf_message(message)
	4105	ISLTPK(I,J)=1
	4106	ENDIF
	4107
	4108
	4109	! cold start: determine liquid soil water content (SH2O)
	4110	! SH2O <= SMC for T < 273.149K (-0.001C)
	4111
	4112	IF (STC(I,K,J) .LT. 273.149) THEN
	4113
	4114	! first guess following explicit solution for Flerchinger Eqn from Koren
	4115	! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
	4116
	4117	BX = BETA(ISLTPK(I,J))
	4118	IF ( BETA(ISLTPK(I,J)) .GT. BLIM ) BX = BLIM
	4119
	4120	if ( GRAV*(-PSIS(ISLTPK(I,J))) .eq. 0 ) then
	4121	write(message,*) 'TROUBLE'
	4122	CALL wrf_message(message)
	4123	write(message,*) 'I,J: ', i,J
	4124	CALL wrf_message(message)
	4125	write(message,*) 'grav, isltpk, psis(isltpk): ', grav,isltpk(I,J),&
	4126	psis(isltpk(I,J))
	4127	CALL wrf_message(message)
	4128	endif
	4129
	4130	if (BX .eq. 0 .or. STC(I,K,J) .eq. 0) then
	4131	write(message,*) 'TROUBLE -- I,J,BX, STC: ', I,J,BX,STC(I,K,J)
	4132	CALL wrf_message(message)
	4133	endif
	4134	FK = (((HLICE/(GRAV(-PSIS(ISLTPK(I,J))))) &
	4135	((STC(I,K,J)-T0)/STC(I,K,J)))** &
	4136	(-1/BX))*SMCMAX(ISLTPK(I,J))
	4137	IF (FK .LT. 0.02) FK = 0.02
	4138	SH2O(I,K,J) = MIN ( FK, SMC(I,K,J) )
	4139	! ----------------------------------------------------------------------
	4140	! now use iterative solution for liquid soil water content using
	4141	! FUNCTION FRH2O (from the Eta "NOAH" land-surface model) with the
	4142	! initial guess for SH2O from above explicit first guess.
	4143
	4144	SH2O(I,K,J)=FRH2O_init(STC(I,K,J),SMC(I,K,J),SH2O(I,K,J), &
	4145	SMCMAX(ISLTPK(I,J)),BETA(ISLTPK(I,J)), &
	4146	PSIS(ISLTPK(I,J)))
	4147
	4148	ELSE ! above freezing
	4149	SH2O(I,K,J)=SMC(I,K,J)
	4150	ENDIF
	4151
	4152
	4153	ELSE ! water point
	4154	SH2O(I,K,J)=SMC(I,K,J)
	4155
	4156	ENDIF ! test on land/ice/sea
	4157	if (SH2O(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
	4158	write(message,*) 'SH2O > THAN SMCMAX ', I,J,SH2O(I,K,J),SMCMAX(ISLTPK(I,J)),SMC(I,K,J)
	4159	CALL wrf_message(message)
	4160	endif
	4161
	4162	ENDDO
	4163	ENDDO
	4164	ENDDO
	4165
	4166	END SUBROUTINE NMM_SH2O
	4167
	4168	!-------------------------------------------------------------------
	4169
	4170	FUNCTION FRH2O_init(TKELV,SMC,SH2O,SMCMAX,B,PSIS)
	4171
	4172	IMPLICIT NONE
	4173
	4174	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4175	! PURPOSE: CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT
	4176	! IF TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION
	4177	! TO SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF
	4178	! KOREN ET AL. (1999, JGR, VOL 104(D16), 19569-19585).
	4179	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4180	!
	4181	! New version (JUNE 2001): much faster and more accurate newton iteration
	4182	! achieved by first taking log of eqn cited above -- less than 4
	4183	! (typically 1 or 2) iterations achieves convergence. Also, explicit
	4184	! 1-step solution option for special case of parameter Ck=0, which reduces
	4185	! the original implicit equation to a simpler explicit form, known as the
	4186	! ""Flerchinger Eqn". Improved handling of solution in the limit of
	4187	! freezing point temperature T0.
	4188	!
	4189	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4190	!
	4191	! INPUT:
	4192	!
	4193	! TKELV.........Temperature (Kelvin)
	4194	! SMC...........Total soil moisture content (volumetric)
	4195	! SH2O..........Liquid soil moisture content (volumetric)
	4196	! SMCMAX........Saturation soil moisture content (from REDPRM)
	4197	! B.............Soil type "B" parameter (from REDPRM)
	4198	! PSIS..........Saturated soil matric potential (from REDPRM)
	4199	!
	4200	! OUTPUT:
	4201	! FRH2O.........supercooled liquid water content.
	4202	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4203
	4204	REAL B
	4205	REAL BLIM
	4206	REAL BX
	4207	REAL CK
	4208	REAL DENOM
	4209	REAL DF
	4210	REAL DH2O
	4211	REAL DICE
	4212	REAL DSWL
	4213	REAL ERROR
	4214	REAL FK
	4215	REAL FRH2O_init
	4216	REAL GS
	4217	REAL HLICE
	4218	REAL PSIS
	4219	REAL SH2O
	4220	REAL SMC
	4221	REAL SMCMAX
	4222	REAL SWL
	4223	REAL SWLK
	4224	REAL TKELV
	4225	REAL T0
	4226
	4227	INTEGER NLOG
	4228	INTEGER KCOUNT
	4229	PARAMETER (CK=8.0)
	4230	! PARAMETER (CK=0.0)
	4231	PARAMETER (BLIM=5.5)
	4232	! PARAMETER (BLIM=7.0)
	4233	PARAMETER (ERROR=0.005)
	4234
	4235	PARAMETER (HLICE=3.335E5)
	4236	PARAMETER (GS = 9.81)
	4237	PARAMETER (DICE=920.0)
	4238	PARAMETER (DH2O=1000.0)
	4239	PARAMETER (T0=273.15)
	4240
	4241	! ### LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM) ####
	4242	! ### SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT ####
	4243	! ### IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES ####
	4244	! ################################################################
	4245	!
	4246	BX = B
	4247	IF ( B .GT. BLIM ) BX = BLIM
	4248	! ------------------------------------------------------------------
	4249
	4250	! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
	4251	NLOG=0
	4252	KCOUNT=0
	4253
	4254	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4255	! C IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
	4256	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4257
	4258
	4259	IF (TKELV .GT. (T0 - 1.E-3)) THEN
	4260
	4261	FRH2O_init=SMC
	4262
	4263	ELSE
	4264
	4265	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4266	IF (CK .NE. 0.0) THEN
	4267
	4268	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4269	! CCCCCCCCC OPTION 1: ITERATED SOLUTION FOR NONZERO CK CCCCCCCCCCC
	4270	! CCCCCCCCCCCC IN KOREN ET AL, JGR, 1999, EQN 17 CCCCCCCCCCCCCCCCC
	4271
	4272	! INITIAL GUESS FOR SWL (frozen content)
	4273	SWL = SMC-SH2O
	4274	! KEEP WITHIN BOUNDS.
	4275	IF (SWL .GT. (SMC-0.02)) SWL=SMC-0.02
	4276	IF(SWL .LT. 0.) SWL=0.
	4277	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4278	! C START OF ITERATIONS
	4279	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4280	DO WHILE (NLOG .LT. 10 .AND. KCOUNT .EQ. 0)
	4281	NLOG = NLOG+1
	4282	DF = ALOG(( PSISGS/HLICE ) ( ( 1.+CKSWL )2. ) &
	4283	( SMCMAX/(SMC-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
	4284	DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( SMC - SWL )
	4285	SWLK = SWL - DF/DENOM
	4286	! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
	4287	IF (SWLK .GT. (SMC-0.02)) SWLK = SMC - 0.02
	4288	IF(SWLK .LT. 0.) SWLK = 0.
	4289	! MATHEMATICAL SOLUTION BOUNDS APPLIED.
	4290	DSWL=ABS(SWLK-SWL)
	4291	SWL=SWLK
	4292	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4293	! CC IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
	4294	! CC WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
	4295	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4296	IF ( DSWL .LE. ERROR ) THEN
	4297	KCOUNT=KCOUNT+1
	4298	END IF
	4299	END DO
	4300	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4301	! C END OF ITERATIONS
	4302	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4303	! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
	4304	FRH2O_init = SMC - SWL
	4305
	4306	! CCCCCCCCCCCCCCCCCCCCCCCC END OPTION 1 CCCCCCCCCCCCCCCCCCCCCCCCCCC
	4307
	4308	ENDIF
	4309
	4310	IF (KCOUNT .EQ. 0) THEN
	4311	! Print*,'Flerchinger used in NEW version. Iterations=',NLOG
	4312
	4313	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	4314	! CCCCC OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0 CCCCCCCC
	4315	! CCCCCCCCCCCCC IN KOREN ET AL., JGR, 1999, EQN 17 CCCCCCCCCCCCCCC
	4316
	4317	FK=(((HLICE/(GS(-PSIS)))((TKELV-T0)/TKELV))*(-1/BX))SMCMAX
	4318	! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
	4319	IF (FK .LT. 0.02) FK = 0.02
	4320	FRH2O_init = MIN ( FK, SMC )
	4321
	4322	! CCCCCCCCCCCCCCCCCCCCCCCCC END OPTION 2 CCCCCCCCCCCCCCCCCCCCCCCCCC
	4323
	4324	ENDIF
	4325
	4326	ENDIF
	4327
	4328	RETURN
	4329
	4330	END FUNCTION FRH2O_init
	4331
	4332
	4333	!--------------------------------------------------------------------
	4334
	4335	SUBROUTINE init_module_initialize
	4336	END SUBROUTINE init_module_initialize
	4337
	4338	!---------------------------------------------------------------------
	4339
	4340	END MODULE module_initialize_real

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

Download in other formats:

Original Format