Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_initialize_real.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 175.3 KB

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

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

Download in other formats:

Original Format