Context Navigation

module_initialize_real.F

Last change on this file 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: 240.6 KB

Line
1	!REAL:MODEL_LAYER:INITIALIZATION
2
3	#ifndef VERT_UNIT
4	! This MODULE holds the routines which are used to perform various initializations
5	! for the individual domains, specifically for the Eulerian, mass-based coordinate.
6
7	!-----------------------------------------------------------------------
8
9	MODULE module_initialize_real
10
11	USE module_bc
12	USE module_configure
13	USE module_domain
14	USE module_io_domain
15	USE module_model_constants
16	USE module_state_description
17	USE module_timing
18	USE module_soil_pre
19	USE module_date_time
20	USE module_llxy
21	#ifdef DM_PARALLEL
22	USE module_dm
23	USE module_comm_dm, ONLY : &
24	HALO_EM_INIT_1_sub &
25	,HALO_EM_INIT_2_sub &
26	,HALO_EM_INIT_3_sub &
27	,HALO_EM_INIT_4_sub &
28	,HALO_EM_INIT_5_sub &
29	,HALO_EM_VINTERP_UV_1_sub
30	#endif
31
32	REAL , SAVE :: p_top_save
33	INTEGER :: internal_time_loop
34
35	CONTAINS
36
37	!-------------------------------------------------------------------
38
39	SUBROUTINE init_domain ( grid )
40
41	IMPLICIT NONE
42
43	! Input space and data. No gridded meteorological data has been stored, though.
44
45	! TYPE (domain), POINTER :: grid
46	TYPE (domain) :: grid
47
48	! Local data.
49
50	INTEGER :: idum1, idum2
51
52	CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
53
54	CALL init_domain_rk( grid &
55	!
56	#include "actual_new_args.inc"
57	!
58	)
59	END SUBROUTINE init_domain
60
61	!-------------------------------------------------------------------
62
63	SUBROUTINE init_domain_rk ( grid &
64	!
65	#include "dummy_new_args.inc"
66	!
67	)
68
69	USE module_optional_input
70	IMPLICIT NONE
71
72	! Input space and data. No gridded meteorological data has been stored, though.
73
74	! TYPE (domain), POINTER :: grid
75	TYPE (domain) :: grid
76
77	#include "dummy_new_decl.inc"
78
79	TYPE (grid_config_rec_type) :: config_flags
80
81	! Local domain indices and counters.
82
83	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
84	INTEGER :: loop , num_seaice_changes
85
86	INTEGER :: ids, ide, jds, jde, kds, kde, &
87	ims, ime, jms, jme, kms, kme, &
88	its, ite, jts, jte, kts, kte, &
89	ips, ipe, jps, jpe, kps, kpe, &
90	i, j, k
91
92	INTEGER :: imsx, imex, jmsx, jmex, kmsx, kmex, &
93	ipsx, ipex, jpsx, jpex, kpsx, kpex, &
94	imsy, imey, jmsy, jmey, kmsy, kmey, &
95	ipsy, ipey, jpsy, jpey, kpsy, kpey
96
97	INTEGER :: ns
98
99	! Local data
100
101	INTEGER :: error
102	INTEGER :: im, num_3d_m, num_3d_s
103	REAL :: p_surf, p_level
104	REAL :: cof1, cof2
105	REAL :: qvf , qvf1 , qvf2 , pd_surf
106	REAL :: p00 , t00 , a , tiso
107	REAL :: hold_znw , ptemp
108	REAL :: vap_pres_mb , sat_vap_pres_mb
109	LOGICAL :: were_bad
110
111	LOGICAL :: stretch_grid, dry_sounding, debug
112	INTEGER IICOUNT
113
114	REAL :: p_top_requested , temp
115	INTEGER :: num_metgrid_levels
116	REAL , DIMENSION(max_eta) :: eta_levels
117	REAL :: max_dz
118
119	! INTEGER , PARAMETER :: nl_max = 1000
120	! REAL , DIMENSION(nl_max) :: grid%dn
121
122	integer::oops1,oops2
123
124	REAL :: zap_close_levels
125	INTEGER :: force_sfc_in_vinterp
126	INTEGER :: interp_type , lagrange_order , extrap_type , t_extrap_type
127	LOGICAL :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
128	LOGICAL :: we_have_tavgsfc , we_have_tsk
129
130	INTEGER :: lev500 , loop_count
131	REAL :: zl , zu , pl , pu , z500 , dz500 , tvsfc , dpmu
132
133	LOGICAL , PARAMETER :: want_full_levels = .TRUE.
134	LOGICAL , PARAMETER :: want_half_levels = .FALSE.
135
136	CHARACTER (LEN=80) :: a_message
137	REAL :: max_mf
138
139	! Excluded middle.
140
141	LOGICAL :: any_valid_points
142	INTEGER :: i_valid , j_valid
143
144	!-- Carsel and Parrish [1988]
145	REAL , DIMENSION(100) :: lqmi
146
147	REAL :: t_start , t_end
148
149	! Dimension information stored in grid data structure.
150
151	CALL cpu_time(t_start)
152	CALL get_ijk_from_grid ( grid , &
153	ids, ide, jds, jde, kds, kde, &
154	ims, ime, jms, jme, kms, kme, &
155	ips, ipe, jps, jpe, kps, kpe, &
156	imsx, imex, jmsx, jmex, kmsx, kmex, &
157	ipsx, ipex, jpsx, jpex, kpsx, kpex, &
158	imsy, imey, jmsy, jmey, kmsy, kmey, &
159	ipsy, ipey, jpsy, jpey, kpsy, kpey )
160	its = ips ; ite = ipe ; jts = jps ; jte = jpe ; kts = kps ; kte = kpe
161
162
163	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
164
165	! Send out a quick message about the time steps based on the map scale factors.
166
167	IF ( ( internal_time_loop .EQ. 1 ) .AND. ( grid%id .EQ. 1 ) .AND. ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) ) THEN
168	max_mf = grid%msft(its,jts)
169	DO j=jts,MIN(jde-1,jte)
170	DO i=its,MIN(ide-1,ite)
171	max_mf = MAX ( max_mf , grid%msft(i,j) )
172	END DO
173	END DO
174	#if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
175	max_mf = wrf_dm_max_real ( max_mf )
176	#endif
177	WRITE ( a_message , FMT='(A,F5.2,A)' ) 'Max map factor in domain 1 = ',max_mf,'. Scale the dt in the model accordingly.'
178	CALL wrf_message ( a_message )
179	END IF
180
181	! Check to see if the boundary conditions are set properly in the namelist file.
182	! This checks for sufficiency and redundancy.
183
184	CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
185
186	! Some sort of "this is the first time" initialization. Who knows.
187
188	grid%step_number = 0
189	grid%itimestep=0
190
191	! Pull in the info in the namelist to compare it to the input data.
192
193	grid%real_data_init_type = model_config_rec%real_data_init_type
194
195	! To define the base state, we call a USER MODIFIED routine to set the three
196	! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K),
197	! and A (temperature difference, from 1000 mb to 300 mb, K).
198
199	CALL const_module_initialize ( p00 , t00 , a , tiso )
200
201	! Save these constants to write out in model output file
202
203	grid%t00 = t00
204	grid%p00 = p00
205	grid%tlp = a
206	grid%tiso = tiso
207
208	! Are there any hold-ups to us bypassing the middle of the domain? These
209	! holdups would be situations where we need data in the middle of the domain.
210	! FOr example, if this si the first time period, we need the full domain
211	! processed for ICs. Also, if there is some sort of gridded FDDA turned on, or
212	! if the SST update is activated, then we can't just blow off the middle of the
213	! domain all willy-nilly. Other cases of these hold-ups? Sure - what if the
214	! user wants to smooth the CG topo, we need several rows and columns available.
215	! What if the lat/lon proj is used, then we need to run a spectral filter on
216	! the topo. Both are killers when trying to ignore data in the middle of the
217	! domain.
218
219	! If hold_ups = .F., then there are no hold-ups to excluding the middle
220	! domain processing. If hold_ups = .T., then there are hold-ups, and we
221	! must process the middle of the domain.
222
223	hold_ups = ( internal_time_loop .EQ. 1 ) .OR. &
224	( config_flags%grid_fdda .NE. 0 ) .OR. &
225	( config_flags%sst_update .EQ. 1 ) .OR. &
226	( config_flags%all_ic_times ) .OR. &
227	( config_flags%smooth_cg_topo ) .OR. &
228	( config_flags%map_proj .EQ. PROJ_CASSINI )
229
230	! There are a few checks that we need to do when the input data comes in with the middle
231	! excluded by WPS.
232
233	IF ( flag_excluded_middle .NE. 0 ) THEN
234
235	! If this time period of data from WPS has the middle excluded, it had better be OK for
236	! us to have a hole.
237
238	IF ( hold_ups ) THEN
239	WRITE ( a_message,* ) 'None of the following are allowed to be TRUE : '
240	CALL wrf_message ( a_message )
241	WRITE ( a_message,* ) ' ( internal_time_loop .EQ. 1 ) ', ( internal_time_loop .EQ. 1 )
242	CALL wrf_message ( a_message )
243	WRITE ( a_message,* ) ' ( config_flags%grid_fdda .NE. 0 ) ', ( config_flags%grid_fdda .NE. 0 )
244	CALL wrf_message ( a_message )
245	WRITE ( a_message,* ) ' ( config_flags%sst_update .EQ. 1 ) ', ( config_flags%sst_update .EQ. 1 )
246	CALL wrf_message ( a_message )
247	WRITE ( a_message,* ) ' ( config_flags%all_ic_times ) ', ( config_flags%all_ic_times )
248	CALL wrf_message ( a_message )
249	WRITE ( a_message,* ) ' ( config_flags%smooth_cg_topo ) ', ( config_flags%smooth_cg_topo )
250	CALL wrf_message ( a_message )
251	WRITE ( a_message,* ) ' ( config_flags%map_proj .EQ. PROJ_CASSINI ) ', ( config_flags%map_proj .EQ. PROJ_CASSINI )
252	CALL wrf_message ( a_message )
253
254	WRITE ( a_message,* ) 'Problems, we cannot have excluded middle data from WPS'
255	CALL wrf_error_fatal ( a_message )
256	END IF
257
258	! Make sure that the excluded middle data from metgrid is "wide enough". We only have to check
259	! when the excluded middle was actually used in WPS.
260
261	IF ( config_flags%spec_bdy_width .GT. flag_excluded_middle ) THEN
262	WRITE ( a_message,* ) 'The WRF &bdy_control namelist.input spec_bdy_width = ', config_flags%spec_bdy_width
263	CALL wrf_message ( a_message )
264	WRITE ( a_message,* ) 'The WPS &metgrid namelist.wps process_only_bdy width = ',flag_excluded_middle
265	CALL wrf_message ( a_message )
266	WRITE ( a_message,* ) 'WPS process_only_bdy must be >= WRF spec_bdy_width'
267	CALL wrf_error_fatal ( a_message )
268	END IF
269	END IF
270	em_width = config_flags%spec_bdy_width
271
272	! We need to find if there are any valid non-excluded-middle points in this
273	! tile. If so, then we need to hang on to a valid i,j location.
274
275	any_valid_points = .false.
276	find_valid : DO j = jts,jte
277	DO i = its,ite
278	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
279	any_valid_points = .true.
280	i_valid = i
281	j_valid = j
282	EXIT find_valid
283	END DO
284	END DO find_valid
285
286	! Replace traditional seaice field with optional seaice (AFWA source)
287
288	IF ( flag_icefrac .EQ. 1 ) THEN
289	DO j=jts,MIN(jde-1,jte)
290	DO i=its,MIN(ide-1,ite)
291	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
292	grid%xice(i,j) = grid%icefrac_gc(i,j)
293	END DO
294	END DO
295	END IF
296
297	! Fix the snow (water equivalent depth, kg/m^2) and the snowh (physical snow
298	! depth, m) fields.
299
300	IF ( ( flag_snow .EQ. 0 ) .AND. ( flag_snowh .EQ. 0 ) ) THEN
301	DO j=jts,MIN(jde-1,jte)
302	DO i=its,MIN(ide-1,ite)
303	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
304	grid%snow(i,j) = 0.
305	grid%snowh(i,j) = 0.
306	END DO
307	END DO
308
309	ELSE IF ( ( flag_snow .EQ. 0 ) .AND. ( flag_snowh .EQ. 1 ) ) THEN
310	DO j=jts,MIN(jde-1,jte)
311	DO i=its,MIN(ide-1,ite)
312	! ( m -> kg/m^2 ) & ( reduce to liquid, 5:1 ratio )
313	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
314	grid%snow(i,j) = grid%snowh(i,j) * 1000. / 5.
315	END DO
316	END DO
317
318	ELSE IF ( ( flag_snow .EQ. 1 ) .AND. ( flag_snowh .EQ. 0 ) ) THEN
319	DO j=jts,MIN(jde-1,jte)
320	DO i=its,MIN(ide-1,ite)
321	! ( kg/m^2 -> m) & ( liquid to snow depth, 5:1 ratio )
322	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
323	grid%snowh(i,j) = grid%snow(i,j) / 1000. * 5.
324	END DO
325	END DO
326
327	END IF
328
329	! For backward compatibility, we might need to assign the map factors from
330	! what they were, to what they are.
331
332	IF ( ( config_flags%polar ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
333	DO j=max(jds+1,jts),min(jde-1,jte)
334	DO i=its,min(ide-1,ite)
335	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
336	grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
337	END DO
338	END DO
339	IF(jts == jds) THEN
340	DO i=its,ite
341	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
342	grid%msfvx(i,jts) = 0.
343	grid%msfvx_inv(i,jts) = 0.
344	END DO
345	END IF
346	IF(jte == jde) THEN
347	DO i=its,ite
348	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
349	grid%msfvx(i,jte) = 0.
350	grid%msfvx_inv(i,jte) = 0.
351	END DO
352	END IF
353	ELSE IF ( ( config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
354	DO j=jts,jte
355	DO i=its,min(ide-1,ite)
356	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
357	grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
358	END DO
359	END DO
360	ELSE IF ( ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .NE. 1 ) ) THEN
361	DO j=jts,jte
362	DO i=its,ite
363	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
364	grid%msfvx(i,j) = grid%msfv(i,j)
365	grid%msfvy(i,j) = grid%msfv(i,j)
366	grid%msfux(i,j) = grid%msfu(i,j)
367	grid%msfuy(i,j) = grid%msfu(i,j)
368	grid%msftx(i,j) = grid%msft(i,j)
369	grid%msfty(i,j) = grid%msft(i,j)
370	ENDDO
371	ENDDO
372	DO j=jts,min(jde,jte)
373	DO i=its,min(ide-1,ite)
374	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
375	grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
376	END DO
377	END DO
378	ELSE IF ( ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
379	IF ( grid%msfvx(its,jts) .EQ. 0 ) THEN
380	CALL wrf_error_fatal ( 'Maybe you do not have the new map factors, try re-running geogrid' )
381	END IF
382	DO j=jts,min(jde,jte)
383	DO i=its,min(ide-1,ite)
384	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
385	grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
386	END DO
387	END DO
388	ELSE IF ( ( config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .NE. 1 ) ) THEN
389	CALL wrf_error_fatal ( 'Neither SI data nor older metgrid data can initialize a global domain' )
390	ENDIF
391
392	! Check to see what available surface temperatures we have.
393
394	IF ( flag_tavgsfc .EQ. 1 ) THEN
395	we_have_tavgsfc = .TRUE.
396	ELSE
397	we_have_tavgsfc = .FALSE.
398	END IF
399
400	IF ( flag_tsk .EQ. 1 ) THEN
401	we_have_tsk = .TRUE.
402	ELSE
403	we_have_tsk = .FALSE.
404	END IF
405
406	IF ( config_flags%use_tavg_for_tsk ) THEN
407	IF ( we_have_tsk .OR. we_have_tavgsfc ) THEN
408	! we are OK
409	ELSE
410	CALL wrf_error_fatal ( 'We either need TSK or TAVGSFC, verify these fields are coming from WPS' )
411	END IF
412
413	! Since we require a skin temperature in the model, we can use the average 2-m temperature if provided.
414
415	IF ( we_have_tavgsfc ) THEN
416	DO j=jts,min(jde,jte)
417	DO i=its,min(ide-1,ite)
418	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
419	grid%tsk(i,j) = grid%tavgsfc(i,j)
420	END DO
421	END DO
422	END IF
423	END IF
424
425	! Is there any vertical interpolation to do? The "old" data comes in on the correct
426	! vertical locations already.
427
428	IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
429
430	! If this is data from the PINTERP program, it is emulating METGRID output.
431	! One of the caveats of this data is the way that the vertical structure is
432	! handled. We take the k=1 level and toss it (it is disposable), and we
433	! swap in the surface data. This is done for all of the 3d fields about
434	! which we show some interest: u, v, t, rh, ght, and p. For u, v, and rh,
435	! we assume no interesting vertical structure, and just assign the 1000 mb
436	! data. We directly use the 2-m temp for surface temp. We use the surface
437	! pressure field and the topography elevation for the lowest level of
438	! pressure and height, respectively.
439
440	IF ( flag_pinterp .EQ. 1 ) THEN
441
442	WRITE ( a_message , * ) 'Data from P_INTERP program, filling k=1 level with artificial surface fields.'
443	CALL wrf_message ( a_message )
444	DO j=jts,jte
445	DO i=its,ite
446	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
447	grid%u_gc(i,1,j) = grid%u_gc(i,2,j)
448	grid%v_gc(i,1,j) = grid%v_gc(i,2,j)
449	grid%rh_gc(i,1,j) = grid%rh_gc(i,2,j)
450	grid%t_gc(i,1,j) = grid%t2(i,j)
451	grid%ght_gc(i,1,j) = grid%ht(i,j)
452	grid%p_gc(i,1,j) = grid%psfc(i,j)
453	END DO
454	END DO
455	flag_psfc = 0
456
457	END IF
458
459	! Variables that are named differently between SI and WPS.
460
461	DO j = jts, MIN(jte,jde-1)
462	DO i = its, MIN(ite,ide-1)
463	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
464	grid%tsk(i,j) = grid%tsk_gc(i,j)
465	grid%tmn(i,j) = grid%tmn_gc(i,j)
466	grid%xlat(i,j) = grid%xlat_gc(i,j)
467	grid%xlong(i,j) = grid%xlong_gc(i,j)
468	grid%ht(i,j) = grid%ht_gc(i,j)
469	END DO
470	END DO
471
472	! A user could request that the most coarse grid has the
473	! topography along the outer boundary smoothed. This smoothing
474	! is similar to the coarse/nest interface. The outer rows and
475	! cols come from the existing large scale topo, and then the
476	! next several rows/cols are a linear ramp of the large scale
477	! model and the hi-res topo from WPS. We only do this for the
478	! coarse grid since we are going to make the interface consistent
479	! in the model betwixt the CG and FG domains.
480
481	IF ( ( config_flags%smooth_cg_topo ) .AND. &
482	( grid%id .EQ. 1 ) .AND. &
483	( flag_soilhgt .EQ. 1) ) THEN
484	CALL blend_terrain ( grid%toposoil , grid%ht , &
485	ids , ide , jds , jde , 1 , 1 , &
486	ims , ime , jms , jme , 1 , 1 , &
487	ips , ipe , jps , jpe , 1 , 1 )
488
489	END IF
490
491	! Filter the input topography if this is a polar projection.
492
493	IF ( ( config_flags%polar ) .AND. ( grid%fft_filter_lat .GT. 90 ) ) THEN
494	CALL wrf_error_fatal ( 'If the polar boundary condition is used, then fft_filter_lat must be set in namelist.input' )
495	END IF
496
497	IF ( config_flags%map_proj .EQ. PROJ_CASSINI ) THEN
498	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
499
500	! We stick the topo and map fac in an unused 3d array. The map scale
501	! factor and computational latitude are passed along for the ride
502	! (part of the transpose process - we only do 3d arrays) to determine
503	! "how many" values are used to compute the mean. We want a number
504	! that is consistent with the original grid resolution.
505
506
507	DO j = jts, MIN(jte,jde-1)
508	DO k = kts, kte
509	DO i = its, MIN(ite,ide-1)
510	grid%t_init(i,k,j) = 1.
511	END DO
512	END DO
513	DO i = its, MIN(ite,ide-1)
514	grid%t_init(i,1,j) = grid%ht(i,j)
515	grid%t_init(i,2,j) = grid%msftx(i,j)
516	grid%t_init(i,3,j) = grid%clat(i,j)
517	END DO
518	END DO
519
520	# include "XPOSE_POLAR_FILTER_TOPO_z2x.inc"
521
522	! Retrieve the 2d arrays for topo, map factors, and the
523	! computational latitude.
524
525	DO j = jpsx, MIN(jpex,jde-1)
526	DO i = ipsx, MIN(ipex,ide-1)
527	grid%ht_xxx(i,j) = grid%t_xxx(i,1,j)
528	grid%mf_xxx(i,j) = grid%t_xxx(i,2,j)
529	grid%clat_xxx(i,j) = grid%t_xxx(i,3,j)
530	END DO
531	END DO
532
533	! Get a mean topo field that is consistent with the grid
534	! distance on each computational latitude loop.
535
536	CALL filter_topo ( grid%ht_xxx , grid%clat_xxx , grid%mf_xxx , &
537	grid%fft_filter_lat , &
538	ids, ide, jds, jde, 1 , 1 , &
539	imsx, imex, jmsx, jmex, 1, 1, &
540	ipsx, ipex, jpsx, jpex, 1, 1 )
541
542	! Stick the filtered topo back into the dummy 3d array to
543	! transpose it back to "all z on a patch".
544
545	DO j = jpsx, MIN(jpex,jde-1)
546	DO i = ipsx, MIN(ipex,ide-1)
547	grid%t_xxx(i,1,j) = grid%ht_xxx(i,j)
548	END DO
549	END DO
550
551	# include "XPOSE_POLAR_FILTER_TOPO_x2z.inc"
552
553	! Get the un-transposed topo data.
554
555	DO j = jts, MIN(jte,jde-1)
556	DO i = its, MIN(ite,ide-1)
557	grid%ht(i,j) = grid%t_init(i,1,j)
558	END DO
559	END DO
560	#else
561	CALL filter_topo ( grid%ht , grid%clat , grid%msftx , &
562	grid%fft_filter_lat , &
563	ids, ide, jds, jde, 1,1, &
564	ims, ime, jms, jme, 1,1, &
565	its, ite, jts, jte, 1,1 )
566	#endif
567	END IF
568
569	! If we have any input low-res surface pressure, we store it.
570
571	IF ( flag_psfc .EQ. 1 ) THEN
572	DO j = jts, MIN(jte,jde-1)
573	DO i = its, MIN(ite,ide-1)
574	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
575	grid%psfc_gc(i,j) = grid%psfc(i,j)
576	grid%p_gc(i,1,j) = grid%psfc(i,j)
577	END DO
578	END DO
579	END IF
580
581	! If we have the low-resolution surface elevation, stick that in the
582	! "input" locations of the 3d height. We still have the "hi-res" topo
583	! stuck in the grid%ht array. The grid%landmask if test is required as some sources
584	! have ZERO elevation over water (thank you very much).
585
586	IF ( flag_soilhgt .EQ. 1) THEN
587	DO j = jts, MIN(jte,jde-1)
588	DO i = its, MIN(ite,ide-1)
589	! IF ( grid%landmask(i,j) .GT. 0.5 ) THEN
590	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
591	grid%ght_gc(i,1,j) = grid%toposoil(i,j)
592	grid%ht_gc(i,j)= grid%toposoil(i,j)
593	! END IF
594	END DO
595	END DO
596	END IF
597
598	! The number of vertical levels in the input data. There is no staggering for
599	! different variables.
600
601	num_metgrid_levels = grid%num_metgrid_levels
602
603	! For UM data, swap incoming extra (theta-based) pressure with the standardly
604	! named (rho-based) pressure.
605
606	IF ( flag_ptheta .EQ. 1 ) THEN
607	DO j = jts, MIN(jte,jde-1)
608	DO k = 1 , num_metgrid_levels
609	DO i = its, MIN(ite,ide-1)
610	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
611	ptemp = grid%p_gc(i,k,j)
612	grid%p_gc(i,k,j) = grid%prho_gc(i,k,j)
613	grid%prho_gc(i,k,j) = ptemp
614	END DO
615	END DO
616	END DO
617
618	! For UM data, the "surface" and the "first hybrid" level for the theta-level data fields are the same.
619	! Average the surface (k=1) and the second hybrid level (k=num_metgrid_levels-1) to get the first hybrid
620	! layer. We only do this for the theta-level data: pressure, temperature, specific humidity, and
621	! geopotential height (i.e. we do not modify u, v, or the rho-based pressure).
622
623	DO j = jts, MIN(jte,jde-1)
624	DO i = its, MIN(ite,ide-1)
625	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
626	grid% p_gc(i,num_metgrid_levels,j) = ( grid% p_gc(i,1,j) + grid% p_gc(i,num_metgrid_levels-1,j) ) * 0.5
627	grid% t_gc(i,num_metgrid_levels,j) = ( grid% t_gc(i,1,j) + grid% t_gc(i,num_metgrid_levels-1,j) ) * 0.5
628	grid% sh_gc(i,num_metgrid_levels,j) = ( grid% sh_gc(i,1,j) + grid% sh_gc(i,num_metgrid_levels-1,j) ) * 0.5
629	grid%ght_gc(i,num_metgrid_levels,j) = ( grid%ght_gc(i,1,j) + grid%ght_gc(i,num_metgrid_levels-1,j) ) * 0.5
630	END DO
631	END DO
632	END IF
633
634	! Check for and semi-fix missing surface fields.
635
636	IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
637	k = 2
638	ELSE
639	k = num_metgrid_levels
640	END IF
641
642	IF ( grid%t_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
643	DO j = jts, MIN(jte,jde-1)
644	DO i = its, MIN(ite,ide-1)
645	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
646	grid%t_gc(i,1,j) = grid%t_gc(i,k,j)
647	END DO
648	END DO
649	config_flags%use_surface = .FALSE.
650	grid%use_surface = .FALSE.
651	WRITE ( a_message , * ) 'Missing surface temp, replaced with closest level, use_surface set to false.'
652	CALL wrf_message ( a_message )
653	END IF
654
655	IF ( grid%rh_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
656	DO j = jts, MIN(jte,jde-1)
657	DO i = its, MIN(ite,ide-1)
658	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
659	grid%rh_gc(i,1,j) = grid%rh_gc(i,k,j)
660	END DO
661	END DO
662	config_flags%use_surface = .FALSE.
663	grid%use_surface = .FALSE.
664	WRITE ( a_message , * ) 'Missing surface RH, replaced with closest level, use_surface set to false.'
665	CALL wrf_message ( a_message )
666	END IF
667
668	IF ( grid%u_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
669	DO j = jts, MIN(jte,jde-1)
670	DO i = its, ite
671	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
672	grid%u_gc(i,1,j) = grid%u_gc(i,k,j)
673	END DO
674	END DO
675	config_flags%use_surface = .FALSE.
676	grid%use_surface = .FALSE.
677	WRITE ( a_message , * ) 'Missing surface u wind, replaced with closest level, use_surface set to false.'
678	CALL wrf_message ( a_message )
679	END IF
680
681	IF ( grid%v_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
682	DO j = jts, jte
683	DO i = its, MIN(ite,ide-1)
684	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
685	grid%v_gc(i,1,j) = grid%v_gc(i,k,j)
686	END DO
687	END DO
688	config_flags%use_surface = .FALSE.
689	grid%use_surface = .FALSE.
690	WRITE ( a_message , * ) 'Missing surface v wind, replaced with closest level, use_surface set to false.'
691	CALL wrf_message ( a_message )
692	END IF
693
694
695	! Compute the mixing ratio from the input relative humidity.
696
697	IF ( ( flag_qv .NE. 1 ) .AND. ( flag_sh .NE. 1 ) ) THEN
698	IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
699	k = 2
700	ELSE
701	k = num_metgrid_levels
702	END IF
703
704	IF ( config_flags%rh2qv_method .eq. 1 ) THEN
705	CALL rh_to_mxrat1(grid%rh_gc, grid%t_gc, grid%p_gc, grid%qv_gc , &
706	config_flags%rh2qv_wrt_liquid , &
707	config_flags%qv_max_p_safe , &
708	config_flags%qv_max_flag , config_flags%qv_max_value , &
709	config_flags%qv_min_p_safe , &
710	config_flags%qv_min_flag , config_flags%qv_min_value , &
711	ids , ide , jds , jde , 1 , num_metgrid_levels , &
712	ims , ime , jms , jme , 1 , num_metgrid_levels , &
713	its , ite , jts , jte , 1 , num_metgrid_levels )
714	ELSE IF ( config_flags%rh2qv_method .eq. 2 ) THEN
715	CALL rh_to_mxrat2(grid%rh_gc, grid%t_gc, grid%p_gc, grid%qv_gc , &
716	config_flags%rh2qv_wrt_liquid , &
717	config_flags%qv_max_p_safe , &
718	config_flags%qv_max_flag , config_flags%qv_max_value , &
719	config_flags%qv_min_p_safe , &
720	config_flags%qv_min_flag , config_flags%qv_min_value , &
721	ids , ide , jds , jde , 1 , num_metgrid_levels , &
722	ims , ime , jms , jme , 1 , num_metgrid_levels , &
723	its , ite , jts , jte , 1 , num_metgrid_levels )
724	END IF
725
726
727	ELSE IF ( flag_sh .EQ. 1 ) THEN
728	IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
729	k = 2
730	ELSE
731	k = num_metgrid_levels
732	END IF
733	IF ( grid%sh_gc(i_valid,kts,j_valid) .LT. 1.e-6 ) THEN
734	DO j = jts, MIN(jte,jde-1)
735	DO i = its, MIN(ite,ide-1)
736	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
737	grid%sh_gc(i,1,j) = grid%sh_gc(i,k,j)
738	END DO
739	END DO
740	END IF
741
742	DO j = jts, MIN(jte,jde-1)
743	DO k = 1 , num_metgrid_levels
744	DO i = its, MIN(ite,ide-1)
745	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
746	grid%qv_gc(i,k,j) = grid%sh_gc(i,k,j) /( 1. - grid%sh_gc(i,k,j) )
747	sat_vap_pres_mb = 0.611210.EXP(17.67*(grid%t_gc(i,k,j)-273.15)/(grid%t_gc(i,k,j)-29.65))
748	vap_pres_mb = grid%qv_gc(i,k,j) * grid%p_gc(i,k,j)/100. / (grid%qv_gc(i,k,j) + 0.622 )
749	IF ( sat_vap_pres_mb .GT. 0 ) THEN
750	grid%rh_gc(i,k,j) = ( vap_pres_mb / sat_vap_pres_mb ) * 100.
751	ELSE
752	grid%rh_gc(i,k,j) = 0.
753	END IF
754	END DO
755	END DO
756	END DO
757
758	END IF
759
760	! Some data sets do not provide a 3d geopotential height field.
761
762	IF ( grid%ght_gc(i_valid,grid%num_metgrid_levels/2,j_valid) .LT. 1 ) THEN
763	DO j = jts, MIN(jte,jde-1)
764	DO k = kts+1 , grid%num_metgrid_levels
765	DO i = its, MIN(ite,ide-1)
766	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
767	grid%ght_gc(i,k,j) = grid%ght_gc(i,k-1,j) - &
768	R_d / g * 0.5 * ( grid%t_gc(i,k ,j) * ( 1 + 0.608 * grid%qv_gc(i,k ,j) ) + &
769	grid%t_gc(i,k-1,j) * ( 1 + 0.608 * grid%qv_gc(i,k-1,j) ) ) * &
770	LOG ( grid%p_gc(i,k,j) / grid%p_gc(i,k-1,j) )
771	END DO
772	END DO
773	END DO
774	END IF
775
776	! Assign surface fields with original input values. If this is hybrid data,
777	! the values are not exactly representative. However - this is only for
778	! plotting purposes and such at the 0h of the forecast, so we are not all that
779	! worried.
780
781	DO j = jts, min(jde-1,jte)
782	DO i = its, min(ide,ite)
783	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
784	grid%u10(i,j)=grid%u_gc(i,1,j)
785	END DO
786	END DO
787
788	DO j = jts, min(jde,jte)
789	DO i = its, min(ide-1,ite)
790	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
791	grid%v10(i,j)=grid%v_gc(i,1,j)
792	END DO
793	END DO
794
795	DO j = jts, min(jde-1,jte)
796	DO i = its, min(ide-1,ite)
797	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
798	grid%t2(i,j)=grid%t_gc(i,1,j)
799	END DO
800	END DO
801
802	IF ( flag_qv .EQ. 1 ) THEN
803	DO j = jts, min(jde-1,jte)
804	DO i = its, min(ide-1,ite)
805	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
806	grid%q2(i,j)=grid%qv_gc(i,1,j)
807	END DO
808	END DO
809	END IF
810
811	! The requested ptop for real data cases.
812
813	p_top_requested = grid%p_top_requested
814
815	! Compute the top pressure, grid%p_top. For isobaric data, this is just the
816	! top level. For the generalized vertical coordinate data, we find the
817	! max pressure on the top level. We have to be careful of two things:
818	! 1) the value has to be communicated, 2) the value can not increase
819	! at subsequent times from the initial value.
820
821	IF ( internal_time_loop .EQ. 1 ) THEN
822	CALL find_p_top ( grid%p_gc , grid%p_top , &
823	ids , ide , jds , jde , 1 , num_metgrid_levels , &
824	ims , ime , jms , jme , 1 , num_metgrid_levels , &
825	its , ite , jts , jte , 1 , num_metgrid_levels )
826
827	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
828	grid%p_top = wrf_dm_max_real ( grid%p_top )
829	#endif
830
831	! Compare the requested grid%p_top with the value available from the input data.
832
833	IF ( p_top_requested .LT. grid%p_top ) THEN
834	print *,'p_top_requested = ',p_top_requested
835	print *,'allowable grid%p_top in data = ',grid%p_top
836	CALL wrf_error_fatal ( 'p_top_requested < grid%p_top possible from data' )
837	END IF
838
839	! The grid%p_top valus is the max of what is available from the data and the
840	! requested value. We have already compared <, so grid%p_top is directly set to
841	! the value in the namelist.
842
843	grid%p_top = p_top_requested
844
845	! For subsequent times, we have to remember what the grid%p_top for the first
846	! time was. Why? If we have a generalized vert coordinate, the grid%p_top value
847	! could fluctuate.
848
849	p_top_save = grid%p_top
850
851	ELSE
852	CALL find_p_top ( grid%p_gc , grid%p_top , &
853	ids , ide , jds , jde , 1 , num_metgrid_levels , &
854	ims , ime , jms , jme , 1 , num_metgrid_levels , &
855	its , ite , jts , jte , 1 , num_metgrid_levels )
856
857	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
858	grid%p_top = wrf_dm_max_real ( grid%p_top )
859	#endif
860	IF ( grid%p_top .GT. p_top_save ) THEN
861	print *,'grid%p_top from last time period = ',p_top_save
862	print *,'grid%p_top from this time period = ',grid%p_top
863	CALL wrf_error_fatal ( 'grid%p_top > previous value' )
864	END IF
865	grid%p_top = p_top_save
866	ENDIF
867
868	! Get the monthly values interpolated to the current date for the traditional monthly
869	! fields of green-ness fraction and background albedo.
870
871	CALL monthly_interp_to_date ( grid%greenfrac , current_date , grid%vegfra , &
872	ids , ide , jds , jde , kds , kde , &
873	ims , ime , jms , jme , kms , kme , &
874	its , ite , jts , jte , kts , kte )
875
876	CALL monthly_interp_to_date ( grid%albedo12m , current_date , grid%albbck , &
877	ids , ide , jds , jde , kds , kde , &
878	ims , ime , jms , jme , kms , kme , &
879	its , ite , jts , jte , kts , kte )
880
881	! Get the min/max of each i,j for the monthly green-ness fraction.
882
883	CALL monthly_min_max ( grid%greenfrac , grid%shdmin , grid%shdmax , &
884	ids , ide , jds , jde , kds , kde , &
885	ims , ime , jms , jme , kms , kme , &
886	its , ite , jts , jte , kts , kte )
887
888	! The model expects the green-ness values in percent, not fraction.
889
890	DO j = jts, MIN(jte,jde-1)
891	DO i = its, MIN(ite,ide-1)
892	grid%vegfra(i,j) = grid%vegfra(i,j) * 100.
893	grid%shdmax(i,j) = grid%shdmax(i,j) * 100.
894	grid%shdmin(i,j) = grid%shdmin(i,j) * 100.
895	END DO
896	END DO
897
898	! The model expects the albedo fields as a fraction, not a percent. Set the
899	! water values to 8%.
900
901	DO j = jts, MIN(jte,jde-1)
902	DO i = its, MIN(ite,ide-1)
903	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
904	grid%albbck(i,j) = grid%albbck(i,j) / 100.
905	grid%snoalb(i,j) = grid%snoalb(i,j) / 100.
906	IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
907	grid%albbck(i,j) = 0.08
908	grid%snoalb(i,j) = 0.08
909	END IF
910	END DO
911	END DO
912
913	! If the pressure levels in the middle of the atmosphere are upside down, then
914	! this is hybrid data. Computing the new surface pressure should use sfcprs2.
915
916	IF ( grid%p_gc(i_valid,num_metgrid_levels/2,j_valid) .LT. grid%p_gc(i_valid,num_metgrid_levels/2+1,j_valid) ) THEN
917	config_flags%sfcp_to_sfcp = .TRUE.
918	END IF
919
920	! Two ways to get the surface pressure. 1) If we have the low-res input surface
921	! pressure and the low-res topography, then we can do a simple hydrostatic
922	! relation. 2) Otherwise we compute the surface pressure from the sea-level
923	! pressure.
924	! Note that on output, grid%psfc is now hi-res. The low-res surface pressure and
925	! elevation are grid%psfc_gc and grid%ht_gc (same as grid%ght_gc(k=1)).
926
927	IF ( ( flag_psfc .EQ. 1 ) .AND. &
928	( flag_soilhgt .EQ. 1 ) .AND. &
929	( flag_slp .EQ. 1 ) .AND. &
930	( .NOT. config_flags%sfcp_to_sfcp ) ) THEN
931	CALL sfcprs3(grid%ght_gc, grid%p_gc, grid%ht, &
932	grid%pslv_gc, grid%psfc, &
933	ids , ide , jds , jde , 1 , num_metgrid_levels , &
934	ims , ime , jms , jme , 1 , num_metgrid_levels , &
935	its , ite , jts , jte , 1 , num_metgrid_levels )
936	ELSE IF ( ( flag_psfc .EQ. 1 ) .AND. &
937	( flag_soilhgt .EQ. 1 ) .AND. &
938	( config_flags%sfcp_to_sfcp ) ) THEN
939	CALL sfcprs2(grid%t_gc, grid%qv_gc, grid%ght_gc, grid%psfc_gc, grid%ht, &
940	grid%tavgsfc, grid%p_gc, grid%psfc, we_have_tavgsfc, &
941	ids , ide , jds , jde , 1 , num_metgrid_levels , &
942	ims , ime , jms , jme , 1 , num_metgrid_levels , &
943	its , ite , jts , jte , 1 , num_metgrid_levels )
944	ELSE IF ( flag_slp .EQ. 1 ) THEN
945	CALL sfcprs (grid%t_gc, grid%qv_gc, grid%ght_gc, grid%pslv_gc, grid%ht, &
946	grid%tavgsfc, grid%p_gc, grid%psfc, we_have_tavgsfc, &
947	ids , ide , jds , jde , 1 , num_metgrid_levels , &
948	ims , ime , jms , jme , 1 , num_metgrid_levels , &
949	its , ite , jts , jte , 1 , num_metgrid_levels )
950	ELSE
951	WRITE(a_message,FMT='(3(A,I2),A,L1)') 'ERROR in psfc: flag_psfc = ',flag_psfc, &
952	', flag_soilhgt = ',flag_soilhgt , &
953	', flag_slp = ',flag_slp , &
954	', sfcp_to_sfcp = ',config_flags%sfcp_to_sfcp
955	CALL wrf_message ( a_message )
956	CALL wrf_error_fatal ( 'not enough info for a p sfc computation' )
957	END IF
958
959	! If we have no input surface pressure, we'd better stick something in there.
960
961	IF ( flag_psfc .NE. 1 ) THEN
962	DO j = jts, MIN(jte,jde-1)
963	DO i = its, MIN(ite,ide-1)
964	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
965	grid%psfc_gc(i,j) = grid%psfc(i,j)
966	grid%p_gc(i,1,j) = grid%psfc(i,j)
967	END DO
968	END DO
969	END IF
970
971	! Integrate the mixing ratio to get the vapor pressure.
972
973	CALL integ_moist ( grid%qv_gc , grid%p_gc , grid%pd_gc , grid%t_gc , grid%ght_gc , grid%intq_gc , &
974	ids , ide , jds , jde , 1 , num_metgrid_levels , &
975	ims , ime , jms , jme , 1 , num_metgrid_levels , &
976	its , ite , jts , jte , 1 , num_metgrid_levels )
977
978	! If this is UM data, the same moisture removed from the "theta" level pressure data can
979	! be removed from the "rho" level pressures. This is an approximation. We'll revisit to
980	! see if this is a bad idea.
981
982	IF ( flag_ptheta .EQ. 1 ) THEN
983	DO j = jts, MIN(jte,jde-1)
984	DO k = num_metgrid_levels-1 , 1 , -1
985	DO i = its, MIN(ite,ide-1)
986	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
987	ptemp = ((grid%p_gc(i,k,j) - grid%pd_gc(i,k,j)) + (grid%p_gc(i,k+1,j) - grid%pd_gc(i,k+1,j)))/2
988	grid%pdrho_gc(i,k,j) = grid%prho_gc(i,k,j) - ptemp
989	END DO
990	END DO
991	END DO
992	END IF
993
994
995	! Compute the difference between the dry, total surface pressure (input) and the
996	! dry top pressure (constant).
997
998	CALL p_dts ( grid%mu0 , grid%intq_gc , grid%psfc , grid%p_top , &
999	ids , ide , jds , jde , 1 , num_metgrid_levels , &
1000	ims , ime , jms , jme , 1 , num_metgrid_levels , &
1001	its , ite , jts , jte , 1 , num_metgrid_levels )
1002
1003	! Compute the dry, hydrostatic surface pressure.
1004
1005	CALL p_dhs ( grid%pdhs , grid%ht , p00 , t00 , a , &
1006	ids , ide , jds , jde , kds , kde , &
1007	ims , ime , jms , jme , kms , kme , &
1008	its , ite , jts , jte , kts , kte )
1009
1010	! Compute the eta levels if not defined already.
1011
1012	IF ( grid%znw(1) .NE. 1.0 ) THEN
1013
1014	eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
1015	max_dz = model_config_rec%max_dz
1016
1017	CALL compute_eta ( grid%znw , &
1018	eta_levels , max_eta , max_dz , &
1019	grid%p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso , &
1020	ids , ide , jds , jde , kds , kde , &
1021	ims , ime , jms , jme , kms , kme , &
1022	its , ite , jts , jte , kts , kte )
1023	END IF
1024
1025	! The input field is temperature, we want potential temp.
1026
1027	CALL t_to_theta ( grid%t_gc , grid%p_gc , p00 , &
1028	ids , ide , jds , jde , 1 , num_metgrid_levels , &
1029	ims , ime , jms , jme , 1 , num_metgrid_levels , &
1030	its , ite , jts , jte , 1 , num_metgrid_levels )
1031
1032	IF ( flag_slp .EQ. 1 ) THEN
1033
1034	! On the eta surfaces, compute the dry pressure = mu eta, stored in
1035	! grid%pb, since it is a pressure, and we don't need another kms:kme 3d
1036	! array floating around. The grid%pb array is re-computed as the base pressure
1037	! later after the vertical interpolations are complete.
1038
1039	CALL p_dry ( grid%mu0 , grid%znw , grid%p_top , grid%pb , want_full_levels , &
1040	ids , ide , jds , jde , kds , kde , &
1041	ims , ime , jms , jme , kms , kme , &
1042	its , ite , jts , jte , kts , kte )
1043
1044	! All of the vertical interpolations are done in dry-pressure space. The
1045	! input data has had the moisture removed (grid%pd_gc). The target levels (grid%pb)
1046	! had the vapor pressure removed from the surface pressure, then they were
1047	! scaled by the eta levels.
1048
1049	interp_type = 2
1050	lagrange_order = grid%lagrange_order
1051	lowest_lev_from_sfc = .FALSE.
1052	use_levels_below_ground = .TRUE.
1053	use_surface = .TRUE.
1054	zap_close_levels = grid%zap_close_levels
1055	force_sfc_in_vinterp = 0
1056	t_extrap_type = grid%t_extrap_type
1057	extrap_type = 1
1058
1059	! For the height field, the lowest level pressure is the slp (approximately "dry"). The
1060	! lowest level of the input height field (to be associated with slp) then is an array
1061	! of zeros.
1062
1063	DO j = jts, MIN(jte,jde-1)
1064	DO i = its, MIN(ite,ide-1)
1065	grid%psfc_gc(i,j) = grid%pd_gc(i,1,j)
1066	grid%pd_gc(i,1,j) = grid%pslv_gc(i,j) - ( grid%p_gc(i,1,j) - grid%pd_gc(i,1,j) )
1067	grid%ht_gc(i,j) = grid%ght_gc(i,1,j)
1068	grid%ght_gc(i,1,j) = 0.
1069	END DO
1070	END DO
1071
1072	CALL vert_interp ( grid%ght_gc , grid%pd_gc , grid%ph0 , grid%pb , &
1073	num_metgrid_levels , 'Z' , &
1074	interp_type , lagrange_order , extrap_type , &
1075	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1076	zap_close_levels , force_sfc_in_vinterp , &
1077	ids , ide , jds , jde , kds , kde , &
1078	ims , ime , jms , jme , kms , kme , &
1079	its , ite , jts , jte , kts , kte )
1080
1081	! Put things back to normal.
1082
1083	DO j = jts, MIN(jte,jde-1)
1084	DO i = its, MIN(ite,ide-1)
1085	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1086	grid%pd_gc(i,1,j) = grid%psfc_gc(i,j)
1087	grid%ght_gc(i,1,j) = grid%ht_gc(i,j)
1088	END DO
1089	END DO
1090
1091	END IF
1092
1093	! Now the rest of the variables on half-levels to inteprolate.
1094
1095	CALL p_dry ( grid%mu0 , grid%znw , grid%p_top , grid%pb , want_half_levels , &
1096	ids , ide , jds , jde , kds , kde , &
1097	ims , ime , jms , jme , kms , kme , &
1098	its , ite , jts , jte , kts , kte )
1099
1100	interp_type = grid%interp_type
1101	lagrange_order = grid%lagrange_order
1102	lowest_lev_from_sfc = grid%lowest_lev_from_sfc
1103	use_levels_below_ground = grid%use_levels_below_ground
1104	use_surface = grid%use_surface
1105	zap_close_levels = grid%zap_close_levels
1106	force_sfc_in_vinterp = grid%force_sfc_in_vinterp
1107	t_extrap_type = grid%t_extrap_type
1108	extrap_type = grid%extrap_type
1109
1110	! Interpolate RH, diagnose Qv later when have temp and pressure. Temporarily
1111	! store this in the u_1 space, for later diagnosis into Qv and stored into moist.
1112
1113	CALL vert_interp ( grid%rh_gc , grid%pd_gc , grid%u_1 , grid%pb , &
1114	num_metgrid_levels , 'Q' , &
1115	interp_type , lagrange_order , extrap_type , &
1116	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1117	zap_close_levels , force_sfc_in_vinterp , &
1118	ids , ide , jds , jde , kds , kde , &
1119	ims , ime , jms , jme , kms , kme , &
1120	its , ite , jts , jte , kts , kte )
1121
1122	CALL vert_interp ( grid%t_gc , grid%pd_gc , grid%t_2 , grid%pb , &
1123	num_metgrid_levels , 'T' , &
1124	interp_type , lagrange_order , t_extrap_type , &
1125	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1126	zap_close_levels , force_sfc_in_vinterp , &
1127	ids , ide , jds , jde , kds , kde , &
1128	ims , ime , jms , jme , kms , kme , &
1129	its , ite , jts , jte , kts , kte )
1130
1131	CALL vert_interp ( grid%p_gc , grid%pd_gc , grid%p , grid%pb , &
1132	num_metgrid_levels , 'T' , &
1133	interp_type , lagrange_order , t_extrap_type , &
1134	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1135	zap_close_levels , force_sfc_in_vinterp , &
1136	ids , ide , jds , jde , kds , kde , &
1137	ims , ime , jms , jme , kms , kme , &
1138	its , ite , jts , jte , kts , kte )
1139
1140	! Do not have full pressure on eta levels, get a first guess at Qv by using
1141	! dry pressure. The use of u_1 (rh) and v_1 (temperature) is temporary.
1142	! We fix the approximation to Qv after the total pressure is available on
1143	! eta surfaces.
1144
1145	grid%v_1 = grid%t_2
1146	CALL theta_to_t ( grid%v_1 , grid%p , p00 , &
1147	ids , ide , jds , jde , kds , kde , &
1148	ims , ime , jms , jme , kms , kme , &
1149	its , ite , jts , jte , kts , kte )
1150
1151	IF ( config_flags%rh2qv_method .eq. 1 ) THEN
1152	CALL rh_to_mxrat1(grid%u_1, grid%v_1, grid%p , moist(:,:,:,P_QV) , &
1153	config_flags%rh2qv_wrt_liquid , &
1154	config_flags%qv_max_p_safe , &
1155	config_flags%qv_max_flag , config_flags%qv_max_value , &
1156	config_flags%qv_min_p_safe , &
1157	config_flags%qv_min_flag , config_flags%qv_min_value , &
1158	ids , ide , jds , jde , kds , kde , &
1159	ims , ime , jms , jme , kms , kme , &
1160	its , ite , jts , jte , kts , kte-1 )
1161	ELSE IF ( config_flags%rh2qv_method .eq. 2 ) THEN
1162	CALL rh_to_mxrat2(grid%u_1, grid%v_1, grid%p , moist(:,:,:,P_QV) , &
1163	config_flags%rh2qv_wrt_liquid , &
1164	config_flags%qv_max_p_safe , &
1165	config_flags%qv_max_flag , config_flags%qv_max_value , &
1166	config_flags%qv_min_p_safe , &
1167	config_flags%qv_min_flag , config_flags%qv_min_value , &
1168	ids , ide , jds , jde , kds , kde , &
1169	ims , ime , jms , jme , kms , kme , &
1170	its , ite , jts , jte , kts , kte-1 )
1171	END IF
1172
1173	num_3d_m = num_moist
1174	num_3d_s = num_scalar
1175
1176	IF ( flag_qr .EQ. 1 ) THEN
1177	DO im = PARAM_FIRST_SCALAR, num_3d_m
1178	IF ( im .EQ. P_QR ) THEN
1179	CALL vert_interp ( grid%qr_gc , grid%pd_gc , moist(:,:,:,P_QR) , grid%pb , &
1180	num_metgrid_levels , 'Q' , &
1181	interp_type , lagrange_order , extrap_type , &
1182	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1183	zap_close_levels , force_sfc_in_vinterp , &
1184	ids , ide , jds , jde , kds , kde , &
1185	ims , ime , jms , jme , kms , kme , &
1186	its , ite , jts , jte , kts , kte )
1187	END IF
1188	END DO
1189	END IF
1190
1191	IF ( flag_qc .EQ. 1 ) THEN
1192	DO im = PARAM_FIRST_SCALAR, num_3d_m
1193	IF ( im .EQ. P_QC ) THEN
1194	CALL vert_interp ( grid%qc_gc , grid%pd_gc , moist(:,:,:,P_QC) , grid%pb , &
1195	num_metgrid_levels , 'Q' , &
1196	interp_type , lagrange_order , extrap_type , &
1197	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1198	zap_close_levels , force_sfc_in_vinterp , &
1199	ids , ide , jds , jde , kds , kde , &
1200	ims , ime , jms , jme , kms , kme , &
1201	its , ite , jts , jte , kts , kte )
1202	END IF
1203	END DO
1204	END IF
1205
1206	IF ( flag_qi .EQ. 1 ) THEN
1207	DO im = PARAM_FIRST_SCALAR, num_3d_m
1208	IF ( im .EQ. P_QI ) THEN
1209	CALL vert_interp ( grid%qi_gc , grid%pd_gc , moist(:,:,:,P_QI) , grid%pb , &
1210	num_metgrid_levels , 'Q' , &
1211	interp_type , lagrange_order , extrap_type , &
1212	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1213	zap_close_levels , force_sfc_in_vinterp , &
1214	ids , ide , jds , jde , kds , kde , &
1215	ims , ime , jms , jme , kms , kme , &
1216	its , ite , jts , jte , kts , kte )
1217	END IF
1218	END DO
1219	END IF
1220
1221	IF ( flag_qs .EQ. 1 ) THEN
1222	DO im = PARAM_FIRST_SCALAR, num_3d_m
1223	IF ( im .EQ. P_QS ) THEN
1224	CALL vert_interp ( grid%qs_gc , grid%pd_gc , moist(:,:,:,P_QS) , grid%pb , &
1225	num_metgrid_levels , 'Q' , &
1226	interp_type , lagrange_order , extrap_type , &
1227	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1228	zap_close_levels , force_sfc_in_vinterp , &
1229	ids , ide , jds , jde , kds , kde , &
1230	ims , ime , jms , jme , kms , kme , &
1231	its , ite , jts , jte , kts , kte )
1232	END IF
1233	END DO
1234	END IF
1235
1236	IF ( flag_qg .EQ. 1 ) THEN
1237	DO im = PARAM_FIRST_SCALAR, num_3d_m
1238	IF ( im .EQ. P_QG ) THEN
1239	CALL vert_interp ( grid%qg_gc , grid%pd_gc , moist(:,:,:,P_QG) , grid%pb , &
1240	num_metgrid_levels , 'Q' , &
1241	interp_type , lagrange_order , extrap_type , &
1242	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1243	zap_close_levels , force_sfc_in_vinterp , &
1244	ids , ide , jds , jde , kds , kde , &
1245	ims , ime , jms , jme , kms , kme , &
1246	its , ite , jts , jte , kts , kte )
1247	END IF
1248	END DO
1249	END IF
1250
1251	IF ( flag_qh .EQ. 1 ) THEN
1252	DO im = PARAM_FIRST_SCALAR, num_3d_m
1253	IF ( im .EQ. P_QH ) THEN
1254	CALL vert_interp ( grid%qh_gc , grid%pd_gc , moist(:,:,:,P_QH) , grid%pb , &
1255	num_metgrid_levels , 'Q' , &
1256	interp_type , lagrange_order , extrap_type , &
1257	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1258	zap_close_levels , force_sfc_in_vinterp , &
1259	ids , ide , jds , jde , kds , kde , &
1260	ims , ime , jms , jme , kms , kme , &
1261	its , ite , jts , jte , kts , kte )
1262	END IF
1263	END DO
1264	END IF
1265
1266	IF ( flag_qni .EQ. 1 ) THEN
1267	DO im = PARAM_FIRST_SCALAR, num_3d_s
1268	IF ( im .EQ. P_QNI ) THEN
1269	CALL vert_interp ( grid%qni_gc , grid%pd_gc , scalar(:,:,:,P_QNI) , grid%pb , &
1270	num_metgrid_levels , 'Q' , &
1271	interp_type , lagrange_order , extrap_type , &
1272	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1273	zap_close_levels , force_sfc_in_vinterp , &
1274	ids , ide , jds , jde , kds , kde , &
1275	ims , ime , jms , jme , kms , kme , &
1276	its , ite , jts , jte , kts , kte )
1277	END IF
1278	END DO
1279	END IF
1280
1281	! If this is UM data, put the dry rho-based pressure back into the dry pressure array.
1282	! Since the dry pressure is no longer needed, no biggy.
1283
1284	IF ( flag_ptheta .EQ. 1 ) THEN
1285	DO j = jts, MIN(jte,jde-1)
1286	DO k = 1 , num_metgrid_levels
1287	DO i = its, MIN(ite,ide-1)
1288	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1289	grid%pd_gc(i,k,j) = grid%prho_gc(i,k,j)
1290	END DO
1291	END DO
1292	END DO
1293	END IF
1294
1295	#ifdef DM_PARALLEL
1296	ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
1297
1298	! For the U and V vertical interpolation, we need the pressure defined
1299	! at both the locations for the horizontal momentum, which we get by
1300	! averaging two pressure values (i and i-1 for U, j and j-1 for V). The
1301	! pressure field on input (grid%pd_gc) and the pressure of the new coordinate
1302	! (grid%pb) are both communicated with an 8 stencil.
1303
1304	# include "HALO_EM_VINTERP_UV_1.inc"
1305	#endif
1306
1307	CALL vert_interp ( grid%u_gc , grid%pd_gc , grid%u_2 , grid%pb , &
1308	num_metgrid_levels , 'U' , &
1309	interp_type , lagrange_order , extrap_type , &
1310	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1311	zap_close_levels , force_sfc_in_vinterp , &
1312	ids , ide , jds , jde , kds , kde , &
1313	ims , ime , jms , jme , kms , kme , &
1314	its , ite , jts , jte , kts , kte )
1315
1316	CALL vert_interp ( grid%v_gc , grid%pd_gc , grid%v_2 , grid%pb , &
1317	num_metgrid_levels , 'V' , &
1318	interp_type , lagrange_order , extrap_type , &
1319	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
1320	zap_close_levels , force_sfc_in_vinterp , &
1321	ids , ide , jds , jde , kds , kde , &
1322	ims , ime , jms , jme , kms , kme , &
1323	its , ite , jts , jte , kts , kte )
1324
1325	END IF ! <----- END OF VERTICAL INTERPOLATION PART ---->
1326
1327	! Set the temperature of the inland lakes to tavgsfc if the temperature is available
1328	! and islake is > num_veg_cat
1329
1330	num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1331	CALL nl_get_iswater ( grid%id , grid%iswater )
1332	CALL nl_get_islake ( grid%id , grid%islake )
1333
1334	IF ( grid%islake < 0 ) THEN
1335	CALL wrf_debug ( 0 , 'Old data, no inland lake information')
1336	ELSE
1337	IF ( we_have_tavgsfc ) THEN
1338
1339	CALL wrf_debug ( 0 , 'Using inland lakes with average surface temperature')
1340	DO j=jts,MIN(jde-1,jte)
1341	DO i=its,MIN(ide-1,ite)
1342	IF ( grid%landusef(i,grid%islake,j) >= 0.5 ) THEN
1343	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1344	grid%sst(i,j) = grid%tavgsfc(i,j)
1345	grid%tsk(i,j) = grid%tavgsfc(i,j)
1346	END IF
1347	END DO
1348	END DO
1349
1350	ELSE ! We don't have tavgsfc
1351
1352	CALL wrf_debug ( 0 , 'No average surface temperature for use with inland lakes')
1353
1354	END IF
1355	DO j=jts,MIN(jde-1,jte)
1356	DO i=its,MIN(ide-1,ite)
1357	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1358	grid%landusef(i,grid%iswater,j) = grid%landusef(i,grid%iswater,j) + &
1359	grid%landusef(i,grid%islake,j)
1360	grid%landusef(i,grid%islake,j) = 0.
1361	END DO
1362	END DO
1363
1364	END IF
1365
1366	! Save the grid%tsk field for later use in the sea ice surface temperature
1367	! for the Noah LSM scheme.
1368
1369	DO j = jts, MIN(jte,jde-1)
1370	DO i = its, MIN(ite,ide-1)
1371	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1372	grid%tsk_save(i,j) = grid%tsk(i,j)
1373	END DO
1374	END DO
1375
1376	! Protect against bad grid%tsk values over water by supplying grid%sst (if it is
1377	! available, and if the grid%sst is reasonable).
1378
1379	DO j = jts, MIN(jde-1,jte)
1380	DO i = its, MIN(ide-1,ite)
1381	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1382	IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
1383	( grid%sst(i,j) .GT. 170. ) .AND. ( grid%sst(i,j) .LT. 400. ) ) THEN
1384	grid%tsk(i,j) = grid%sst(i,j)
1385	ENDIF
1386	END DO
1387	END DO
1388
1389	! Take the data from the input file and store it in the variables that
1390	! use the WRF naming and ordering conventions.
1391
1392	DO j = jts, MIN(jte,jde-1)
1393	DO i = its, MIN(ite,ide-1)
1394	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1395	IF ( grid%snow(i,j) .GE. 10. ) then
1396	grid%snowc(i,j) = 1.
1397	ELSE
1398	grid%snowc(i,j) = 0.0
1399	END IF
1400	END DO
1401	END DO
1402
1403	! Set flag integers for presence of snowh and soilw fields
1404
1405	grid%ifndsnowh = flag_snowh
1406	IF (num_sw_levels_input .GE. 1) THEN
1407	grid%ifndsoilw = 1
1408	ELSE
1409	grid%ifndsoilw = 0
1410	END IF
1411
1412	! We require input data for the various LSM schemes.
1413
1414	enough_data : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1415
1416	CASE (LSMSCHEME)
1417	IF ( num_st_levels_input .LT. 2 ) THEN
1418	CALL wrf_error_fatal ( 'Not enough soil temperature data for Noah LSM scheme.')
1419	END IF
1420
1421	CASE (RUCLSMSCHEME)
1422	IF ( num_st_levels_input .LT. 2 ) THEN
1423	CALL wrf_error_fatal ( 'Not enough soil temperature data for RUC LSM scheme.')
1424	END IF
1425
1426	CASE (PXLSMSCHEME)
1427	IF ( num_st_levels_input .LT. 2 ) THEN
1428	CALL wrf_error_fatal ( 'Not enough soil temperature data for P-X LSM scheme.')
1429	END IF
1430
1431	END SELECT enough_data
1432
1433	interpolate_soil_tmw : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1434
1435	CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME, PXLSMSCHEME )
1436	CALL process_soil_real ( grid%tsk , grid%tmn , grid%tavgsfc, &
1437	grid%landmask , grid%sst , grid%ht, grid%toposoil, &
1438	st_input , sm_input , sw_input , &
1439	st_levels_input , sm_levels_input , sw_levels_input , &
1440	grid%zs , grid%dzs , grid%tslb , grid%smois , grid%sh2o , &
1441	flag_sst , flag_tavgsfc, &
1442	flag_soilhgt, flag_soil_layers, flag_soil_levels, &
1443	ids , ide , jds , jde , kds , kde , &
1444	ims , ime , jms , jme , kms , kme , &
1445	its , ite , jts , jte , kts , kte , &
1446	model_config_rec%sf_surface_physics(grid%id) , &
1447	model_config_rec%num_soil_layers , &
1448	model_config_rec%real_data_init_type , &
1449	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1450	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
1451
1452	END SELECT interpolate_soil_tmw
1453
1454	! surface_input_source=1 => use data from static file (fractional category as input)
1455	! surface_input_source=2 => use data from grib file (dominant category as input)
1456	! surface_input_source=3 => use dominant data from static file (dominant category as input)
1457
1458	IF ( config_flags%surface_input_source .EQ. 1 ) THEN
1459
1460	! Generate the vegetation and soil category information from the fractional input
1461	! data, or use the existing dominant category fields if they exist.
1462
1463	grid%vegcat (its,jts) = 0
1464	grid%soilcat(its,jts) = 0
1465
1466	num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1467	num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1468	num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1469
1470	CALL process_percent_cat_new ( grid%landmask , &
1471	grid%landusef , grid%soilctop , grid%soilcbot , &
1472	grid%isltyp , grid%ivgtyp , &
1473	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1474	ids , ide , jds , jde , kds , kde , &
1475	ims , ime , jms , jme , kms , kme , &
1476	its , ite , jts , jte , kts , kte , &
1477	model_config_rec%iswater(grid%id) )
1478
1479	! Make all the veg/soil parms the same so as not to confuse the developer.
1480
1481	DO j = jts , MIN(jde-1,jte)
1482	DO i = its , MIN(ide-1,ite)
1483	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1484	grid%vegcat(i,j) = grid%ivgtyp(i,j)
1485	grid%soilcat(i,j) = grid%isltyp(i,j)
1486	END DO
1487	END DO
1488
1489	ELSE IF ( config_flags%surface_input_source .EQ. 2 ) THEN
1490
1491	! Do we have dominant soil and veg data from the input already?
1492
1493	IF ( grid%soilcat(i_valid,j_valid) .GT. 0.5 ) THEN
1494	DO j = jts, MIN(jde-1,jte)
1495	DO i = its, MIN(ide-1,ite)
1496	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1497	grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
1498	END DO
1499	END DO
1500	END IF
1501	IF ( grid%vegcat(i_valid,j_valid) .GT. 0.5 ) THEN
1502	DO j = jts, MIN(jde-1,jte)
1503	DO i = its, MIN(ide-1,ite)
1504	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1505	grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
1506	END DO
1507	END DO
1508	END IF
1509
1510	ELSE IF ( config_flags%surface_input_source .EQ. 3 ) THEN
1511
1512	! Do we have dominant soil and veg data from the static input already?
1513
1514	IF ( grid%sct_dom_gc(i_valid,j_valid) .GT. 0.5 ) THEN
1515	DO j = jts, MIN(jde-1,jte)
1516	DO i = its, MIN(ide-1,ite)
1517	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1518	grid%isltyp(i,j) = NINT( grid%sct_dom_gc(i,j) )
1519	grid%soilcat(i,j) = grid%isltyp(i,j)
1520	END DO
1521	END DO
1522	ELSE
1523	WRITE ( a_message , * ) 'You have set surface_input_source = 3, but your geogrid data does not have valid dominant soil data.'
1524	CALL wrf_error_fatal ( a_message )
1525	END IF
1526	IF ( grid%lu_index(i_valid,j_valid) .GT. 0.5 ) THEN
1527	DO j = jts, MIN(jde-1,jte)
1528	DO i = its, MIN(ide-1,ite)
1529	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1530	grid%ivgtyp(i,j) = NINT( grid%lu_index(i,j) )
1531	grid%vegcat(i,j) = grid%ivgtyp(i,j)
1532	END DO
1533	END DO
1534	ELSE
1535	WRITE ( a_message , * ) 'You have set surface_input_source = 3, but your geogrid data does not have valid dominant land use data.'
1536	CALL wrf_error_fatal ( a_message )
1537	END IF
1538
1539	END IF
1540
1541	! Adjustments for the seaice field PRIOR to the grid%tslb computations. This is
1542	! is for the 5-layer scheme.
1543
1544	num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1545	num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1546	num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1547	CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
1548	CALL nl_get_isice ( grid%id , grid%isice )
1549	CALL nl_get_iswater ( grid%id , grid%iswater )
1550	CALL adjust_for_seaice_pre ( grid%xice , grid%landmask , grid%tsk , grid%ivgtyp , grid%vegcat , grid%lu_index , &
1551	grid%xland , grid%landusef , grid%isltyp , grid%soilcat , grid%soilctop , &
1552	grid%soilcbot , grid%tmn , &
1553	grid%seaice_threshold , &
1554	config_flags%fractional_seaice, &
1555	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1556	grid%iswater , grid%isice , &
1557	model_config_rec%sf_surface_physics(grid%id) , &
1558	ids , ide , jds , jde , kds , kde , &
1559	ims , ime , jms , jme , kms , kme , &
1560	its , ite , jts , jte , kts , kte )
1561
1562	! Land use assignment.
1563
1564	DO j = jts, MIN(jde-1,jte)
1565	DO i = its, MIN(ide-1,ite)
1566	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1567	grid%lu_index(i,j) = grid%ivgtyp(i,j)
1568	IF ( grid%lu_index(i,j) .NE. model_config_rec%iswater(grid%id) ) THEN
1569	grid%landmask(i,j) = 1
1570	grid%xland(i,j) = 1
1571	ELSE
1572	grid%landmask(i,j) = 0
1573	grid%xland(i,j) = 2
1574	END IF
1575	END DO
1576	END DO
1577
1578
1579	! Fix grid%tmn and grid%tsk.
1580
1581	fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1582
1583	CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME, PXLSMSCHEME )
1584	DO j = jts, MIN(jde-1,jte)
1585	DO i = its, MIN(ide-1,ite)
1586	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1587	IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
1588	( grid%sst(i,j) .GT. 170. ) .AND. ( grid%sst(i,j) .LT. 400. ) ) THEN
1589	grid%tmn(i,j) = grid%sst(i,j)
1590	grid%tsk(i,j) = grid%sst(i,j)
1591	ELSE IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
1592	grid%tmn(i,j) = grid%tsk(i,j)
1593	END IF
1594	END DO
1595	END DO
1596	END SELECT fix_tsk_tmn
1597
1598	! Is the grid%tsk reasonable?
1599
1600	IF ( internal_time_loop .NE. 1 ) THEN
1601	DO j = jts, MIN(jde-1,jte)
1602	DO i = its, MIN(ide-1,ite)
1603	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1604	IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
1605	grid%tsk(i,j) = grid%t_2(i,1,j)
1606	END IF
1607	END DO
1608	END DO
1609	ELSE
1610	DO j = jts, MIN(jde-1,jte)
1611	DO i = its, MIN(ide-1,ite)
1612	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1613	IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
1614	print *,'error in the grid%tsk'
1615	print *,'i,j=',i,j
1616	print *,'grid%landmask=',grid%landmask(i,j)
1617	print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1618	if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
1619	grid%tsk(i,j)=grid%tmn(i,j)
1620	else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1621	grid%tsk(i,j)=grid%sst(i,j)
1622	else
1623	CALL wrf_error_fatal ( 'grid%tsk unreasonable' )
1624	end if
1625	END IF
1626	END DO
1627	END DO
1628	END IF
1629
1630	! Is the grid%tmn reasonable?
1631
1632	DO j = jts, MIN(jde-1,jte)
1633	DO i = its, MIN(ide-1,ite)
1634	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1635	IF ( ( ( grid%tmn(i,j) .LT. 170. ) .OR. ( grid%tmn(i,j) .GT. 400. ) ) &
1636	.AND. ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
1637	IF ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) THEN
1638	print *,'error in the grid%tmn'
1639	print *,'i,j=',i,j
1640	print *,'grid%landmask=',grid%landmask(i,j)
1641	print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1642	END IF
1643
1644	if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
1645	grid%tmn(i,j)=grid%tsk(i,j)
1646	else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1647	grid%tmn(i,j)=grid%sst(i,j)
1648	else
1649	CALL wrf_error_fatal ( 'grid%tmn unreasonable' )
1650	endif
1651	END IF
1652	END DO
1653	END DO
1654
1655
1656	! Minimum soil values, residual, from RUC LSM scheme. For input from Noah or EC, and using
1657	! RUC LSM scheme, this must be subtracted from the input total soil moisture. For
1658	! input RUC data and using the Noah LSM scheme, this value must be added to the soil
1659	! moisture input.
1660
1661	lqmi(1:num_soil_top_cat) = &
1662	(/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
1663	0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
1664	0.004, 0.065 /)
1665	! 0.004, 0.065, 0.020, 0.004, 0.008 /) ! has extra levels for playa, lava, and white sand
1666
1667	! At the initial time we care about values of soil moisture and temperature, other times are
1668	! ignored by the model, so we ignore them, too.
1669
1670	IF ( domain_ClockIsStartTime(grid) ) THEN
1671	account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1672
1673	CASE ( LSMSCHEME )
1674	iicount = 0
1675	IF ( flag_soil_layers == 1 ) THEN
1676	DO j = jts, MIN(jde-1,jte)
1677	DO i = its, MIN(ide-1,ite)
1678	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1679	IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 170 ) .and. &
1680	( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
1681	print *,'Noah -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
1682	iicount = iicount + 1
1683	grid%smois(i,:,j) = 0.005
1684	END IF
1685	END DO
1686	END DO
1687	IF ( iicount .GT. 0 ) THEN
1688	print *,'Noah -> Noah: total number of small soil moisture locations = ',iicount
1689	END IF
1690	ELSE IF ( flag_soil_levels == 1 ) THEN
1691	DO j = jts, MIN(jde-1,jte)
1692	DO i = its, MIN(ide-1,ite)
1693	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1694	grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) + lqmi(grid%isltyp(i,j)) , 0.005 )
1695	END DO
1696	END DO
1697	DO j = jts, MIN(jde-1,jte)
1698	DO i = its, MIN(ide-1,ite)
1699	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1700	IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 170 ) .and. &
1701	( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
1702	print *,'RUC -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
1703	iicount = iicount + 1
1704	grid%smois(i,:,j) = 0.005
1705	END IF
1706	END DO
1707	END DO
1708	IF ( iicount .GT. 0 ) THEN
1709	print *,'RUC -> Noah: total number of small soil moisture locations = ',iicount
1710	END IF
1711	END IF
1712
1713	CASE ( RUCLSMSCHEME )
1714	iicount = 0
1715	IF ( flag_soil_layers == 1 ) THEN
1716	DO j = jts, MIN(jde-1,jte)
1717	DO i = its, MIN(ide-1,ite)
1718	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1719	grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) - lqmi(grid%isltyp(i,j)) , 0.005 )
1720	END DO
1721	END DO
1722	ELSE IF ( flag_soil_levels == 1 ) THEN
1723	! no op
1724	END IF
1725
1726	CASE ( PXLSMSCHEME )
1727	iicount = 0
1728	IF ( flag_soil_layers == 1 ) THEN
1729	! no op
1730	ELSE IF ( flag_soil_levels == 1 ) THEN
1731	DO j = jts, MIN(jde-1,jte)
1732	DO i = its, MIN(ide-1,ite)
1733	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1734	grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) + lqmi(grid%isltyp(i,j)) , 0.005 )
1735	END DO
1736	END DO
1737	END IF
1738
1739	END SELECT account_for_zero_soil_moisture
1740	END IF
1741
1742	! Is the grid%tslb reasonable?
1743
1744	IF ( internal_time_loop .NE. 1 ) THEN
1745	DO j = jts, MIN(jde-1,jte)
1746	DO ns = 1 , model_config_rec%num_soil_layers
1747	DO i = its, MIN(ide-1,ite)
1748	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1749	IF ( grid%tslb(i,ns,j) .LT. 170 .or. grid%tslb(i,ns,j) .GT. 400. ) THEN
1750	grid%tslb(i,ns,j) = grid%t_2(i,1,j)
1751	grid%smois(i,ns,j) = 0.3
1752	END IF
1753	END DO
1754	END DO
1755	END DO
1756	ELSE
1757	DO j = jts, MIN(jde-1,jte)
1758	DO i = its, MIN(ide-1,ite)
1759	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1760	IF ( ( ( grid%tslb(i,1,j) .LT. 170. ) .OR. ( grid%tslb(i,1,j) .GT. 400. ) ) .AND. &
1761	( grid%landmask(i,j) .GT. 0.5 ) ) THEN
1762	IF ( ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) .AND. &
1763	( model_config_rec%sf_surface_physics(grid%id) .NE. RUCLSMSCHEME ).AND. &
1764	( model_config_rec%sf_surface_physics(grid%id) .NE. PXLSMSCHEME ) ) THEN
1765	print *,'error in the grid%tslb'
1766	print *,'i,j=',i,j
1767	print *,'grid%landmask=',grid%landmask(i,j)
1768	print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1769	print *,'grid%tslb = ',grid%tslb(i,:,j)
1770	print *,'old grid%smois = ',grid%smois(i,:,j)
1771	grid%smois(i,1,j) = 0.3
1772	grid%smois(i,2,j) = 0.3
1773	grid%smois(i,3,j) = 0.3
1774	grid%smois(i,4,j) = 0.3
1775	END IF
1776
1777	IF ( (grid%tsk(i,j).GT.170. .AND. grid%tsk(i,j).LT.400.) .AND. &
1778	(grid%tmn(i,j).GT.170. .AND. grid%tmn(i,j).LT.400.) ) THEN
1779	fake_soil_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1780	CASE ( SLABSCHEME )
1781	DO ns = 1 , model_config_rec%num_soil_layers
1782	grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
1783	grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
1784	END DO
1785	CASE ( LSMSCHEME , RUCLSMSCHEME, PXLSMSCHEME )
1786	CALL wrf_error_fatal ( 'Assigned constant soil moisture to 0.3, stopping')
1787	DO ns = 1 , model_config_rec%num_soil_layers
1788	grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
1789	grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
1790	END DO
1791	END SELECT fake_soil_temp
1792	else if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
1793	CALL wrf_error_fatal ( 'grid%tslb unreasonable 1' )
1794	DO ns = 1 , model_config_rec%num_soil_layers
1795	grid%tslb(i,ns,j)=grid%tsk(i,j)
1796	END DO
1797	else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1798	CALL wrf_error_fatal ( 'grid%tslb unreasonable 2' )
1799	DO ns = 1 , model_config_rec%num_soil_layers
1800	grid%tslb(i,ns,j)=grid%sst(i,j)
1801	END DO
1802	else if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
1803	CALL wrf_error_fatal ( 'grid%tslb unreasonable 3' )
1804	DO ns = 1 , model_config_rec%num_soil_layers
1805	grid%tslb(i,ns,j)=grid%tmn(i,j)
1806	END DO
1807	else
1808	CALL wrf_error_fatal ( 'grid%tslb unreasonable 4' )
1809	endif
1810	END IF
1811	END DO
1812	END DO
1813	END IF
1814
1815	! Adjustments for the seaice field AFTER the grid%tslb computations. This is
1816	! is for the Noah LSM scheme.
1817
1818	num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1819	num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1820	num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1821	CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
1822	CALL nl_get_isice ( grid%id , grid%isice )
1823	CALL nl_get_iswater ( grid%id , grid%iswater )
1824	CALL adjust_for_seaice_post ( grid%xice , grid%landmask , grid%tsk , grid%tsk_save , &
1825	grid%ivgtyp , grid%vegcat , grid%lu_index , &
1826	grid%xland , grid%landusef , grid%isltyp , grid%soilcat , &
1827	grid%soilctop , &
1828	grid%soilcbot , grid%tmn , grid%vegfra , &
1829	grid%tslb , grid%smois , grid%sh2o , &
1830	grid%seaice_threshold , &
1831	grid%sst,flag_sst, &
1832	config_flags%fractional_seaice, &
1833	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1834	model_config_rec%num_soil_layers , &
1835	grid%iswater , grid%isice , &
1836	model_config_rec%sf_surface_physics(grid%id) , &
1837	ids , ide , jds , jde , kds , kde , &
1838	ims , ime , jms , jme , kms , kme , &
1839	its , ite , jts , jte , kts , kte )
1840
1841	! Let us make sure (again) that the grid%landmask and the veg/soil categories match.
1842
1843	oops1=0
1844	oops2=0
1845	DO j = jts, MIN(jde-1,jte)
1846	DO i = its, MIN(ide-1,ite)
1847	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1848	IF ( ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. &
1849	( grid%ivgtyp(i,j) .NE. config_flags%iswater .OR. grid%isltyp(i,j) .NE. 14 ) ) .OR. &
1850	( ( grid%landmask(i,j) .GT. 0.5 ) .AND. &
1851	( grid%ivgtyp(i,j) .EQ. config_flags%iswater .OR. grid%isltyp(i,j) .EQ. 14 ) ) ) THEN
1852	IF ( grid%tslb(i,1,j) .GT. 1. ) THEN
1853	oops1=oops1+1
1854	grid%ivgtyp(i,j) = 5
1855	grid%isltyp(i,j) = 8
1856	grid%landmask(i,j) = 1
1857	grid%xland(i,j) = 1
1858	ELSE IF ( grid%sst(i,j) .GT. 1. ) THEN
1859	oops2=oops2+1
1860	grid%ivgtyp(i,j) = config_flags%iswater
1861	grid%isltyp(i,j) = 14
1862	grid%landmask(i,j) = 0
1863	grid%xland(i,j) = 2
1864	ELSE
1865	print *,'the grid%landmask and soil/veg cats do not match'
1866	print *,'i,j=',i,j
1867	print *,'grid%landmask=',grid%landmask(i,j)
1868	print *,'grid%ivgtyp=',grid%ivgtyp(i,j)
1869	print *,'grid%isltyp=',grid%isltyp(i,j)
1870	print *,'iswater=', config_flags%iswater
1871	print *,'grid%tslb=',grid%tslb(i,:,j)
1872	print *,'grid%sst=',grid%sst(i,j)
1873	CALL wrf_error_fatal ( 'mismatch_landmask_ivgtyp' )
1874	END IF
1875	END IF
1876	END DO
1877	END DO
1878	if (oops1.gt.0) then
1879	print *,'points artificially set to land : ',oops1
1880	endif
1881	if(oops2.gt.0) then
1882	print *,'points artificially set to water: ',oops2
1883	endif
1884	! fill grid%sst array with grid%tsk if missing in real input (needed for time-varying grid%sst in wrf)
1885	DO j = jts, MIN(jde-1,jte)
1886	DO i = its, MIN(ide-1,ite)
1887	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1888	IF ( flag_sst .NE. 1 ) THEN
1889	grid%sst(i,j) = grid%tsk(i,j)
1890	ENDIF
1891	END DO
1892	END DO
1893	!tgs set snoalb to land value if the water point is covered with ice
1894	DO j = jts, MIN(jde-1,jte)
1895	DO i = its, MIN(ide-1,ite)
1896	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1897	IF ( grid%ivgtyp(i,j) .EQ. config_flags%isice) THEN
1898	grid%snoalb(i,j) = 0.75
1899	ENDIF
1900	END DO
1901	END DO
1902
1903	! From the full level data, we can get the half levels, reciprocals, and layer
1904	! thicknesses. These are all defined at half level locations, so one less level.
1905	! We allow the vertical coordinate to accidently come in upside down. We want
1906	! the first full level to be the ground surface.
1907
1908	! Check whether grid%znw (full level) data are truly full levels. If not, we need to adjust them
1909	! to be full levels.
1910	! in this test, we check if grid%znw(1) is neither 0 nor 1 (within a tolerance of 10**-5)
1911
1912	were_bad = .false.
1913	IF ( ( (grid%znw(1).LT.(1-1.E-5) ) .OR. ( grid%znw(1).GT.(1+1.E-5) ) ).AND. &
1914	( (grid%znw(1).LT.(0-1.E-5) ) .OR. ( grid%znw(1).GT.(0+1.E-5) ) ) ) THEN
1915	were_bad = .true.
1916	print *,'Your grid%znw input values are probably half-levels. '
1917	print *,grid%znw
1918	print *,'WRF expects grid%znw values to be full levels. '
1919	print *,'Adjusting now to full levels...'
1920	! We want to ignore the first value if it's negative
1921	IF (grid%znw(1).LT.0) THEN
1922	grid%znw(1)=0
1923	END IF
1924	DO k=2,kde
1925	grid%znw(k)=2*grid%znw(k)-grid%znw(k-1)
1926	END DO
1927	END IF
1928
1929	! Let's check our changes
1930
1931	IF ( ( ( grid%znw(1) .LT. (1-1.E-5) ) .OR. ( grid%znw(1) .GT. (1+1.E-5) ) ).AND. &
1932	( ( grid%znw(1) .LT. (0-1.E-5) ) .OR. ( grid%znw(1) .GT. (0+1.E-5) ) ) ) THEN
1933	print *,'The input grid%znw height values were half-levels or erroneous. '
1934	print *,'Attempts to treat the values as half-levels and change them '
1935	print *,'to valid full levels failed.'
1936	CALL wrf_error_fatal("bad grid%znw values from input files")
1937	ELSE IF ( were_bad ) THEN
1938	print *,'...adjusted. grid%znw array now contains full eta level values. '
1939	ENDIF
1940
1941	IF ( grid%znw(1) .LT. grid%znw(kde) ) THEN
1942	DO k=1, kde/2
1943	hold_znw = grid%znw(k)
1944	grid%znw(k)=grid%znw(kde+1-k)
1945	grid%znw(kde+1-k)=hold_znw
1946	END DO
1947	END IF
1948
1949	DO k=1, kde-1
1950	grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
1951	grid%rdnw(k) = 1./grid%dnw(k)
1952	grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
1953	END DO
1954
1955	! Now the same sort of computations with the half eta levels, even ANOTHER
1956	! level less than the one above.
1957
1958	DO k=2, kde-1
1959	grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
1960	grid%rdn(k) = 1./grid%dn(k)
1961	grid%fnp(k) = .5* grid%dnw(k )/grid%dn(k)
1962	grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
1963	END DO
1964
1965	! Scads of vertical coefficients.
1966
1967	cof1 = (2.grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))grid%dnw(1)/grid%dn(2)
1968	cof2 = grid%dn(2) /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
1969
1970	grid%cf1 = grid%fnp(2) + cof1
1971	grid%cf2 = grid%fnm(2) - cof1 - cof2
1972	grid%cf3 = cof2
1973
1974	grid%cfn = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
1975	grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
1976
1977	! Inverse grid distances.
1978
1979	grid%rdx = 1./config_flags%dx
1980	grid%rdy = 1./config_flags%dy
1981
1982	! Some of the many weird geopotential initializations that we'll see today: grid%ph0 is total,
1983	! and grid%ph_2 is a perturbation from the base state geopotential. We set the base geopotential
1984	! at the lowest level to terrain elevation * gravity.
1985
1986	DO j=jts,jte
1987	DO i=its,ite
1988	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1989	grid%ph0(i,1,j) = grid%ht(i,j) * g
1990	grid%ph_2(i,1,j) = 0.
1991	END DO
1992	END DO
1993
1994	! Base state potential temperature and inverse density (alpha = 1/rho) from
1995	! the half eta levels and the base-profile surface pressure. Compute 1/rho
1996	! from equation of state. The potential temperature is a perturbation from t0.
1997
1998	DO j = jts, MIN(jte,jde-1)
1999	DO i = its, MIN(ite,ide-1)
2000
2001	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2002
2003	! Base state pressure is a function of eta level and terrain, only, plus
2004	! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
2005	! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
2006
2007	p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
2008
2009
2010	DO k = 1, kte-1
2011	grid%php(i,k,j) = grid%znw(k)*(p_surf - grid%p_top) + grid%p_top ! temporary, full lev base pressure
2012	grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
2013	temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
2014	! temp = t00 + A*LOG(grid%pb(i,k,j)/p00)
2015	grid%t_init(i,k,j) = temp(p00/grid%pb(i,k,j))*(r_d/cp) - t0
2016	grid%alb(i,k,j) = (r_d/p1000mb)(grid%t_init(i,k,j)+t0)(grid%pb(i,k,j)/p1000mb)**cvpm
2017	END DO
2018
2019	! Base state mu is defined as base state surface pressure minus grid%p_top
2020
2021	grid%mub(i,j) = p_surf - grid%p_top
2022
2023	! Dry surface pressure is defined as the following (this mu is from the input file
2024	! computed from the dry pressure). Here the dry pressure is just reconstituted.
2025
2026	pd_surf = grid%mu0(i,j) + grid%p_top
2027
2028	! Integrate base geopotential, starting at terrain elevation. This assures that
2029	! the base state is in exact hydrostatic balance with respect to the model equations.
2030	! This field is on full levels.
2031
2032	grid%phb(i,1,j) = grid%ht(i,j) * g
2033	DO k = 2,kte
2034	grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)grid%mub(i,j)grid%alb(i,k-1,j)
2035	END DO
2036	END DO
2037	END DO
2038
2039	! Fill in the outer rows and columns to allow us to be sloppy.
2040
2041	IF ( ite .EQ. ide ) THEN
2042	i = ide
2043	DO j = jts, MIN(jde-1,jte)
2044	grid%mub(i,j) = grid%mub(i-1,j)
2045	grid%mu_2(i,j) = grid%mu_2(i-1,j)
2046	DO k = 1, kte-1
2047	grid%pb(i,k,j) = grid%pb(i-1,k,j)
2048	grid%t_init(i,k,j) = grid%t_init(i-1,k,j)
2049	grid%alb(i,k,j) = grid%alb(i-1,k,j)
2050	END DO
2051	DO k = 1, kte
2052	grid%phb(i,k,j) = grid%phb(i-1,k,j)
2053	END DO
2054	END DO
2055	END IF
2056
2057	IF ( jte .EQ. jde ) THEN
2058	j = jde
2059	DO i = its, ite
2060	grid%mub(i,j) = grid%mub(i,j-1)
2061	grid%mu_2(i,j) = grid%mu_2(i,j-1)
2062	DO k = 1, kte-1
2063	grid%pb(i,k,j) = grid%pb(i,k,j-1)
2064	grid%t_init(i,k,j) = grid%t_init(i,k,j-1)
2065	grid%alb(i,k,j) = grid%alb(i,k,j-1)
2066	END DO
2067	DO k = 1, kte
2068	grid%phb(i,k,j) = grid%phb(i,k,j-1)
2069	END DO
2070	END DO
2071	END IF
2072
2073	! Compute the perturbation dry pressure (grid%mub + grid%mu_2 + ptop = dry grid%psfc).
2074
2075	DO j = jts, min(jde-1,jte)
2076	DO i = its, min(ide-1,ite)
2077	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2078	grid%mu_2(i,j) = grid%mu0(i,j) - grid%mub(i,j)
2079	END DO
2080	END DO
2081
2082	! Fill in the outer rows and columns to allow us to be sloppy.
2083
2084	IF ( ite .EQ. ide ) THEN
2085	i = ide
2086	DO j = jts, MIN(jde-1,jte)
2087	grid%mu_2(i,j) = grid%mu_2(i-1,j)
2088	END DO
2089	END IF
2090
2091	IF ( jte .EQ. jde ) THEN
2092	j = jde
2093	DO i = its, ite
2094	grid%mu_2(i,j) = grid%mu_2(i,j-1)
2095	END DO
2096	END IF
2097
2098	lev500 = 0
2099	DO j = jts, min(jde-1,jte)
2100	DO i = its, min(ide-1,ite)
2101	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2102
2103	! Assign the potential temperature (perturbation from t0) and qv on all the mass
2104	! point locations.
2105
2106	DO k = 1 , kde-1
2107	grid%t_2(i,k,j) = grid%t_2(i,k,j) - t0
2108	END DO
2109
2110	dpmu = 10001.
2111	loop_count = 0
2112
2113	DO WHILE ( ( ABS(dpmu) .GT. 10. ) .AND. &
2114	( loop_count .LT. 5 ) )
2115
2116	loop_count = loop_count + 1
2117
2118	! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
2119	! equation) down from the top to get the pressure perturbation. First get the pressure
2120	! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
2121
2122	k = kte-1
2123
2124	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
2125	qvf2 = 1./(1.+qvf1)
2126	qvf1 = qvf1*qvf2
2127
2128	grid%p(i,k,j) = - 0.5(grid%mu_2(i,j)+qvf1grid%mub(i,j))/grid%rdnw(k)/qvf2
2129	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
2130	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_2(i,k,j)+t0)qvf&
2131	(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)*cvpm)
2132	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
2133	grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
2134
2135	! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
2136	! inverse density fields (total and perturbation).
2137
2138	DO k=kte-2,1,-1
2139	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
2140	qvf2 = 1./(1.+qvf1)
2141	qvf1 = qvf1*qvf2
2142	grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_2(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1)
2143	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
2144	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_2(i,k,j)+t0)qvf* &
2145	(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
2146	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
2147	grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
2148	END DO
2149
2150	#if 1
2151	! This is the hydrostatic equation used in the model after the small timesteps. In
2152	! the model, grid%al (inverse density) is computed from the geopotential.
2153
2154	DO k = 2,kte
2155	grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
2156	grid%dnw(k-1) * ( (grid%mub(i,j)+grid%mu_2(i,j))*grid%al(i,k-1,j) &
2157	+ grid%mu_2(i,j)*grid%alb(i,k-1,j) )
2158	grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
2159	END DO
2160	#else
2161	! Get the perturbation geopotential from the 3d height array from WPS.
2162
2163	DO k = 2,kte
2164	grid%ph_2(i,k,j) = grid%ph0(i,k,j)*g - grid%phb(i,k,j)
2165	END DO
2166	#endif
2167
2168	! Adjust the column pressure so that the computed 500 mb height is close to the
2169	! input value (of course, not when we are doing hybrid input).
2170
2171	IF ( ( flag_metgrid .EQ. 1 ) .AND. ( i .EQ. i_valid ) .AND. ( j .EQ. j_valid ) ) THEN
2172	DO k = 1 , num_metgrid_levels
2173	IF ( ABS ( grid%p_gc(i,k,j) - 50000. ) .LT. 1. ) THEN
2174	lev500 = k
2175	EXIT
2176	END IF
2177	END DO
2178	END IF
2179
2180	! We only do the adjustment of height if we have the input data on pressure
2181	! surfaces, and folks have asked to do this option.
2182
2183	IF ( ( flag_metgrid .EQ. 1 ) .AND. &
2184	( flag_ptheta .EQ. 1 ) .AND. &
2185	( config_flags%adjust_heights ) .AND. &
2186	( lev500 .NE. 0 ) ) THEN
2187
2188	DO k = 2 , kte-1
2189
2190	! Get the pressures on the full eta levels (grid%php is defined above as
2191	! the full-lev base pressure, an easy array to use for 3d space).
2192
2193	pl = grid%php(i,k ,j) + &
2194	( grid%p(i,k-1 ,j) * ( grid%znw(k ) - grid%znu(k ) ) + &
2195	grid%p(i,k ,j) * ( grid%znu(k-1 ) - grid%znw(k ) ) ) / &
2196	( grid%znu(k-1 ) - grid%znu(k ) )
2197	pu = grid%php(i,k+1,j) + &
2198	( grid%p(i,k-1+1,j) * ( grid%znw(k +1) - grid%znu(k+1) ) + &
2199	grid%p(i,k +1,j) * ( grid%znu(k-1+1) - grid%znw(k+1) ) ) / &
2200	( grid%znu(k-1+1) - grid%znu(k+1) )
2201
2202	! If these pressure levels trap 500 mb, use them to interpolate
2203	! to the 500 mb level of the computed height.
2204
2205	IF ( ( pl .GE. 50000. ) .AND. ( pu .LT. 50000. ) ) THEN
2206	zl = ( grid%ph_2(i,k ,j) + grid%phb(i,k ,j) ) / g
2207	zu = ( grid%ph_2(i,k+1,j) + grid%phb(i,k+1,j) ) / g
2208
2209	z500 = ( zl * ( LOG(50000.) - LOG(pu ) ) + &
2210	zu * ( LOG(pl ) - LOG(50000.) ) ) / &
2211	( LOG(pl) - LOG(pu) )
2212	! z500 = ( zl * ( (50000.) - (pu ) ) + &
2213	! zu * ( (pl ) - (50000.) ) ) / &
2214	! ( (pl) - (pu) )
2215
2216	! Compute the difference of the 500 mb heights (computed minus input), and
2217	! then the change in grid%mu_2. The grid%php is still full-levels, base pressure.
2218
2219	dz500 = z500 - grid%ght_gc(i,lev500,j)
2220	tvsfc = ((grid%t_2(i,1,j)+t0)((grid%p(i,1,j)+grid%php(i,1,j))/p1000mb)(r_d/cp)) &
2221	(1.+0.6*moist(i,1,j,P_QV))
2222	dpmu = ( grid%php(i,1,j) + grid%p(i,1,j) ) * EXP ( g * dz500 / ( r_d * tvsfc ) )
2223	dpmu = dpmu - ( grid%php(i,1,j) + grid%p(i,1,j) )
2224	grid%mu_2(i,j) = grid%mu_2(i,j) - dpmu
2225	EXIT
2226	END IF
2227
2228	END DO
2229	ELSE
2230	dpmu = 0.
2231	END IF
2232
2233	END DO
2234
2235	END DO
2236	END DO
2237
2238	! If this is data from the SI, then we probably do not have the original
2239	! surface data laying around. Note that these are all the lowest levels
2240	! of the respective 3d arrays. For surface pressure, we assume that the
2241	! vertical gradient of grid%p prime is zilch. This is not all that important.
2242	! These are filled in so that the various plotting routines have something
2243	! to play with at the initial time for the model.
2244
2245	IF ( flag_metgrid .NE. 1 ) THEN
2246	DO j = jts, min(jde-1,jte)
2247	DO i = its, min(ide,ite)
2248	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2249	grid%u10(i,j)=grid%u_2(i,1,j)
2250	END DO
2251	END DO
2252
2253	DO j = jts, min(jde,jte)
2254	DO i = its, min(ide-1,ite)
2255	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2256	grid%v10(i,j)=grid%v_2(i,1,j)
2257	END DO
2258	END DO
2259
2260	DO j = jts, min(jde-1,jte)
2261	DO i = its, min(ide-1,ite)
2262	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2263	p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
2264	grid%psfc(i,j)=p_surf + grid%p(i,1,j)
2265	grid%q2(i,j)=moist(i,1,j,P_QV)
2266	grid%th2(i,j)=grid%t_2(i,1,j)+300.
2267	grid%t2(i,j)=grid%th2(i,j)(((grid%p(i,1,j)+grid%pb(i,1,j))/p00)*(r_d/cp))
2268	END DO
2269	END DO
2270
2271	! If this data is from WPS, then we have previously assigned the surface
2272	! data for u, v, and t. If we have an input qv, welp, we assigned that one,
2273	! too. Now we pick up the left overs, and if RH came in - we assign the
2274	! mixing ratio.
2275
2276	ELSE IF ( flag_metgrid .EQ. 1 ) THEN
2277
2278	DO j = jts, min(jde-1,jte)
2279	DO i = its, min(ide-1,ite)
2280	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2281	p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
2282	! grid%psfc(i,j)=p_surf + grid%p(i,1,j)
2283	grid%th2(i,j)=grid%t2(i,j)(p00/(grid%p(i,1,j)+grid%pb(i,1,j)))*(r_d/cp)
2284	END DO
2285	END DO
2286	IF ( flag_qv .NE. 1 ) THEN
2287	DO j = jts, min(jde-1,jte)
2288	DO i = its, min(ide-1,ite)
2289	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2290	grid%q2(i,j)=moist(i,1,j,P_QV)
2291	END DO
2292	END DO
2293	END IF
2294
2295	END IF
2296	CALL cpu_time(t_end)
2297
2298	! Set flag to denote that we are saving original values of HT, MUB, and
2299	! PHB for 2-way nesting and cycling.
2300
2301	grid%save_topo_from_real=1
2302
2303	ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
2304	#ifdef DM_PARALLEL
2305	# include "HALO_EM_INIT_1.inc"
2306	# include "HALO_EM_INIT_2.inc"
2307	# include "HALO_EM_INIT_3.inc"
2308	# include "HALO_EM_INIT_4.inc"
2309	# include "HALO_EM_INIT_5.inc"
2310	#endif
2311
2312	RETURN
2313
2314	END SUBROUTINE init_domain_rk
2315
2316	!---------------------------------------------------------------------
2317
2318	SUBROUTINE const_module_initialize ( p00 , t00 , a , tiso )
2319	USE module_configure
2320	IMPLICIT NONE
2321	! For the real-data-cases only.
2322	REAL , INTENT(OUT) :: p00 , t00 , a , tiso
2323	CALL nl_get_base_pres ( 1 , p00 )
2324	CALL nl_get_base_temp ( 1 , t00 )
2325	CALL nl_get_base_lapse ( 1 , a )
2326	CALL nl_get_iso_temp ( 1 , tiso )
2327	END SUBROUTINE const_module_initialize
2328
2329	!-------------------------------------------------------------------
2330
2331	SUBROUTINE rebalance_driver ( grid )
2332
2333	IMPLICIT NONE
2334
2335	TYPE (domain) :: grid
2336
2337	CALL rebalance( grid &
2338	!
2339	#include "actual_new_args.inc"
2340	!
2341	)
2342
2343	END SUBROUTINE rebalance_driver
2344
2345	!---------------------------------------------------------------------
2346
2347	SUBROUTINE rebalance ( grid &
2348	!
2349	#include "dummy_new_args.inc"
2350	!
2351	)
2352	IMPLICIT NONE
2353
2354	TYPE (domain) :: grid
2355
2356	#include "dummy_new_decl.inc"
2357
2358	TYPE (grid_config_rec_type) :: config_flags
2359
2360	REAL :: p_surf , pd_surf, p_surf_int , pb_int , ht_hold
2361	REAL :: qvf , qvf1 , qvf2
2362	REAL :: p00 , t00 , a , tiso
2363	REAL , DIMENSION(:,:,:) , ALLOCATABLE :: t_init_int
2364
2365	! Local domain indices and counters.
2366
2367	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
2368
2369	INTEGER :: &
2370	ids, ide, jds, jde, kds, kde, &
2371	ims, ime, jms, jme, kms, kme, &
2372	its, ite, jts, jte, kts, kte, &
2373	ips, ipe, jps, jpe, kps, kpe, &
2374	i, j, k
2375
2376	REAL :: temp, temp_int
2377
2378	SELECT CASE ( model_data_order )
2379	CASE ( DATA_ORDER_ZXY )
2380	kds = grid%sd31 ; kde = grid%ed31 ;
2381	ids = grid%sd32 ; ide = grid%ed32 ;
2382	jds = grid%sd33 ; jde = grid%ed33 ;
2383
2384	kms = grid%sm31 ; kme = grid%em31 ;
2385	ims = grid%sm32 ; ime = grid%em32 ;
2386	jms = grid%sm33 ; jme = grid%em33 ;
2387
2388	kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch
2389	its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch
2390	jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
2391
2392	CASE ( DATA_ORDER_XYZ )
2393	ids = grid%sd31 ; ide = grid%ed31 ;
2394	jds = grid%sd32 ; jde = grid%ed32 ;
2395	kds = grid%sd33 ; kde = grid%ed33 ;
2396
2397	ims = grid%sm31 ; ime = grid%em31 ;
2398	jms = grid%sm32 ; jme = grid%em32 ;
2399	kms = grid%sm33 ; kme = grid%em33 ;
2400
2401	its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
2402	jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch
2403	kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch
2404
2405	CASE ( DATA_ORDER_XZY )
2406	ids = grid%sd31 ; ide = grid%ed31 ;
2407	kds = grid%sd32 ; kde = grid%ed32 ;
2408	jds = grid%sd33 ; jde = grid%ed33 ;
2409
2410	ims = grid%sm31 ; ime = grid%em31 ;
2411	kms = grid%sm32 ; kme = grid%em32 ;
2412	jms = grid%sm33 ; jme = grid%em33 ;
2413
2414	its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
2415	kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch
2416	jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
2417
2418	END SELECT
2419
2420	ALLOCATE ( t_init_int(ims:ime,kms:kme,jms:jme) )
2421
2422	! Fill config_flags the options for a particular domain
2423
2424	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
2425
2426	! Some of the many weird geopotential initializations that we'll see today: grid%ph0 is total,
2427	! and grid%ph_2 is a perturbation from the base state geopotential. We set the base geopotential
2428	! at the lowest level to terrain elevation * gravity.
2429
2430	DO j=jts,jte
2431	DO i=its,ite
2432	grid%ph0(i,1,j) = grid%ht_fine(i,j) * g
2433	grid%ph_2(i,1,j) = 0.
2434	END DO
2435	END DO
2436
2437	! To define the base state, we call a USER MODIFIED routine to set the three
2438	! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K),
2439	! and A (temperature difference, from 1000 mb to 300 mb, K), and constant stratosphere
2440	! temp (tiso, K) either from input file or from namelist (for backward compatibiliy).
2441
2442	IF ( config_flags%use_baseparam_fr_nml ) then
2443	! get these from namelist
2444	CALL wrf_message('ndown: using namelist constants')
2445	CALL const_module_initialize ( p00 , t00 , a , tiso )
2446	ELSE
2447	! get these constants from model data
2448	CALL wrf_message('ndown: using constants from file')
2449	t00 = grid%t00
2450	p00 = grid%p00
2451	a = grid%tlp
2452	tiso = grid%tiso
2453
2454	IF (t00 .LT. 100. .or. p00 .LT. 10000.) THEN
2455	WRITE(wrf_err_message,*)&
2456	'ndown_em: did not find base state parameters in wrfout. Add use_baseparam_fr_nml = .t. in &dynamics and rerun'
2457	CALL wrf_error_fatal(TRIM(wrf_err_message))
2458	ENDIF
2459	ENDIF
2460
2461	! Base state potential temperature and inverse density (alpha = 1/rho) from
2462	! the half eta levels and the base-profile surface pressure. Compute 1/rho
2463	! from equation of state. The potential temperature is a perturbation from t0.
2464
2465	DO j = jts, MIN(jte,jde-1)
2466	DO i = its, MIN(ite,ide-1)
2467
2468	! Base state pressure is a function of eta level and terrain, only, plus
2469	! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
2470	! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
2471	! The fine grid terrain is ht_fine, the interpolated is grid%ht.
2472
2473	p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht_fine(i,j)/a/r_d ) *0.5 )
2474	p_surf_int = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j) /a/r_d ) *0.5 )
2475
2476	DO k = 1, kte-1
2477	grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
2478	pb_int = grid%znu(k)*(p_surf_int - grid%p_top) + grid%p_top
2479	temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
2480	! temp = t00 + A*LOG(pb/p00)
2481	grid%t_init(i,k,j) = temp(p00/grid%pb(i,k,j))*(r_d/cp) - t0
2482	! grid%t_init(i,k,j) = (t00 + ALOG(grid%pb(i,k,j)/p00))(p00/grid%pb(i,k,j))**(r_d/cp) - t0
2483	temp_int = MAX ( tiso, t00 + A*LOG(pb_int /p00) )
2484	t_init_int(i,k,j)= temp_int(p00/pb_int )*(r_d/cp) - t0
2485	! t_init_int(i,k,j)= (t00 + ALOG(pb_int /p00))(p00/pb_int )**(r_d/cp) - t0
2486	grid%alb(i,k,j) = (r_d/p1000mb)(grid%t_init(i,k,j)+t0)(grid%pb(i,k,j)/p1000mb)**cvpm
2487	END DO
2488
2489	! Base state mu is defined as base state surface pressure minus grid%p_top
2490
2491	grid%mub(i,j) = p_surf - grid%p_top
2492
2493	! Dry surface pressure is defined as the following (this mu is from the input file
2494	! computed from the dry pressure). Here the dry pressure is just reconstituted.
2495
2496	pd_surf = ( grid%mub(i,j) + grid%mu_2(i,j) ) + grid%p_top
2497
2498	! Integrate base geopotential, starting at terrain elevation. This assures that
2499	! the base state is in exact hydrostatic balance with respect to the model equations.
2500	! This field is on full levels.
2501
2502	grid%phb(i,1,j) = grid%ht_fine(i,j) * g
2503	DO k = 2,kte
2504	grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)grid%mub(i,j)grid%alb(i,k-1,j)
2505	END DO
2506	END DO
2507	END DO
2508
2509	! Replace interpolated terrain with fine grid values.
2510
2511	DO j = jts, MIN(jte,jde-1)
2512	DO i = its, MIN(ite,ide-1)
2513	grid%ht(i,j) = grid%ht_fine(i,j)
2514	END DO
2515	END DO
2516
2517	! Perturbation fields.
2518
2519	DO j = jts, min(jde-1,jte)
2520	DO i = its, min(ide-1,ite)
2521
2522	! The potential temperature is THETAnest = THETAinterp + ( TBARnest - TBARinterp)
2523
2524	DO k = 1 , kde-1
2525	grid%t_2(i,k,j) = grid%t_2(i,k,j) + ( grid%t_init(i,k,j) - t_init_int(i,k,j) )
2526	END DO
2527
2528	! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
2529	! equation) down from the top to get the pressure perturbation. First get the pressure
2530	! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
2531
2532	k = kte-1
2533
2534	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
2535	qvf2 = 1./(1.+qvf1)
2536	qvf1 = qvf1*qvf2
2537
2538	grid%p(i,k,j) = - 0.5(grid%mu_2(i,j)+qvf1grid%mub(i,j))/grid%rdnw(k)/qvf2
2539	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
2540	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_2(i,k,j)+t0)qvf* &
2541	(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
2542	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
2543
2544	! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
2545	! inverse density fields (total and perturbation).
2546
2547	DO k=kte-2,1,-1
2548	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
2549	qvf2 = 1./(1.+qvf1)
2550	qvf1 = qvf1*qvf2
2551	grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_2(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1)
2552	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
2553	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_2(i,k,j)+t0)qvf* &
2554	(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
2555	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
2556	END DO
2557
2558	! This is the hydrostatic equation used in the model after the small timesteps. In
2559	! the model, grid%al (inverse density) is computed from the geopotential.
2560
2561	DO k = 2,kte
2562	grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
2563	grid%dnw(k-1) * ( (grid%mub(i,j)+grid%mu_2(i,j))*grid%al(i,k-1,j) &
2564	+ grid%mu_2(i,j)*grid%alb(i,k-1,j) )
2565	grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
2566	END DO
2567
2568	END DO
2569	END DO
2570
2571	DEALLOCATE ( t_init_int )
2572
2573	ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
2574	#ifdef DM_PARALLEL
2575	# include "HALO_EM_INIT_1.inc"
2576	# include "HALO_EM_INIT_2.inc"
2577	# include "HALO_EM_INIT_3.inc"
2578	# include "HALO_EM_INIT_4.inc"
2579	# include "HALO_EM_INIT_5.inc"
2580	#endif
2581	END SUBROUTINE rebalance
2582
2583	!---------------------------------------------------------------------
2584
2585	RECURSIVE SUBROUTINE find_my_parent ( grid_ptr_in , grid_ptr_out , id_i_am , id_wanted , found_the_id )
2586
2587	! RAR - Modified to correct problem in which the parent of a child domain could
2588	! not be found in the namelist. This condition typically occurs while using the
2589	! "allow_grid" namelist option when an inactive domain comes before an active
2590	! domain in the list, i.e., the domain number of the active domain is greater than
2591	! that of an inactive domain at the same level.
2592	!
2593	USE module_domain
2594
2595	TYPE(domain) , POINTER :: grid_ptr_in , grid_ptr_out
2596	TYPE(domain) , POINTER :: grid_ptr_sibling
2597	INTEGER :: id_wanted , id_i_am
2598	INTEGER :: nest ! RAR
2599	LOGICAL :: found_the_id
2600
2601	found_the_id = .FALSE.
2602	grid_ptr_sibling => grid_ptr_in
2603	nest = 0 ! RAR
2604
2605	DO WHILE ( ASSOCIATED ( grid_ptr_sibling ) )
2606
2607	IF ( grid_ptr_sibling%grid_id .EQ. id_wanted ) THEN
2608	found_the_id = .TRUE.
2609	grid_ptr_out => grid_ptr_sibling
2610	RETURN
2611	! RAR ELSE IF ( grid_ptr_sibling%num_nests .GT. 0 ) THEN
2612	ELSE IF ( grid_ptr_sibling%num_nests .GT. 0 .AND. nest .LT. grid_ptr_sibling%num_nests ) THEN
2613	nest = nest + 1 ! RAR
2614	grid_ptr_sibling => grid_ptr_sibling%nests(nest)%ptr ! RAR
2615	CALL find_my_parent ( grid_ptr_sibling , grid_ptr_out , id_i_am , id_wanted , found_the_id )
2616	IF (.NOT. found_the_id) grid_ptr_sibling => grid_ptr_sibling%parents(1)%ptr ! RAR
2617	ELSE
2618	grid_ptr_sibling => grid_ptr_sibling%sibling
2619	END IF
2620
2621	END DO
2622
2623	END SUBROUTINE find_my_parent
2624
2625	!---------------------------------------------------------------------
2626
2627	RECURSIVE SUBROUTINE find_my_parent2 ( grid_ptr_in , grid_ptr_out , id_wanted , found_the_id )
2628
2629	USE module_domain
2630
2631	TYPE(domain) , POINTER :: grid_ptr_in
2632	TYPE(domain) , POINTER :: grid_ptr_out
2633	INTEGER , INTENT(IN ) :: id_wanted
2634	LOGICAL , INTENT(OUT) :: found_the_id
2635
2636	! Local
2637
2638	TYPE(domain) , POINTER :: grid_ptr_holder
2639	INTEGER :: kid
2640
2641	! Initializations
2642
2643	found_the_id = .FALSE.
2644	grid_ptr_holder => grid_ptr_in
2645
2646
2647	! Have we found the correct location? If so, we can just pop back up with
2648	! the pointer to the right location (i.e. the parent), thank you very much.
2649
2650	IF ( id_wanted .EQ. grid_ptr_in%grid_id ) THEN
2651
2652	found_the_id = .TRUE.
2653	grid_ptr_out => grid_ptr_in
2654
2655
2656	! We gotta keep looking.
2657
2658	ELSE
2659
2660	! We drill down and process each nest from this domain. We don't have to
2661	! worry about siblings, as we are running over all of the kids for this parent,
2662	! so it amounts to the same set of domains being tested.
2663
2664	loop_over_all_kids : DO kid = 1 , grid_ptr_in%num_nests
2665
2666	IF ( ASSOCIATED ( grid_ptr_in%nests(kid)%ptr ) ) THEN
2667
2668	CALL find_my_parent2 ( grid_ptr_in%nests(kid)%ptr , grid_ptr_out , id_wanted , found_the_id )
2669	IF ( found_the_id ) THEN
2670	EXIT loop_over_all_kids
2671	END IF
2672
2673	END IF
2674	END DO loop_over_all_kids
2675
2676	END IF
2677
2678	END SUBROUTINE find_my_parent2
2679
2680	#endif
2681
2682	!---------------------------------------------------------------------
2683
2684	#ifdef VERT_UNIT
2685
2686	!This is a main program for a small unit test for the vertical interpolation.
2687
2688	program vint
2689
2690	implicit none
2691
2692	integer , parameter :: ij = 3
2693	integer , parameter :: keta = 30
2694	integer , parameter :: kgen =20
2695
2696	integer :: ids , ide , jds , jde , kds , kde , &
2697	ims , ime , jms , jme , kms , kme , &
2698	its , ite , jts , jte , kts , kte
2699
2700	integer :: generic
2701
2702	real , dimension(1:ij,kgen,1:ij) :: fo , po
2703	real , dimension(1:ij,1:keta,1:ij) :: fn_calc , fn_interp , pn
2704
2705	integer, parameter :: interp_type = 1 ! 2
2706	! integer, parameter :: lagrange_order = 2 ! 1
2707	integer :: lagrange_order
2708	logical, parameter :: lowest_lev_from_sfc = .FALSE. ! .TRUE.
2709	logical, parameter :: use_levels_below_ground = .FALSE. ! .TRUE.
2710	logical, parameter :: use_surface = .FALSE. ! .TRUE.
2711	real , parameter :: zap_close_levels = 500. ! 100.
2712	integer, parameter :: force_sfc_in_vinterp = 0 ! 6
2713
2714	integer :: k
2715
2716	ids = 1 ; ide = ij ; jds = 1 ; jde = ij ; kds = 1 ; kde = keta
2717	ims = 1 ; ime = ij ; jms = 1 ; jme = ij ; kms = 1 ; kme = keta
2718	its = 1 ; ite = ij ; jts = 1 ; jte = ij ; kts = 1 ; kte = keta
2719
2720	generic = kgen
2721
2722	print *,' '
2723	print *,'------------------------------------'
2724	print *,'UNIT TEST FOR VERTICAL INTERPOLATION'
2725	print *,'------------------------------------'
2726	print *,' '
2727	do lagrange_order = 1 , 2
2728	print *,' '
2729	print *,'------------------------------------'
2730	print *,'Lagrange Order = ',lagrange_order
2731	print *,'------------------------------------'
2732	print *,' '
2733	call fillitup ( fo , po , fn_calc , pn , &
2734	ids , ide , jds , jde , kds , kde , &
2735	ims , ime , jms , jme , kms , kme , &
2736	its , ite , jts , jte , kts , kte , &
2737	generic , lagrange_order )
2738
2739	print *,' '
2740	print *,'Level Pressure Field'
2741	print *,' (Pa) (generic)'
2742	print *,'------------------------------------'
2743	print *,' '
2744	do k = 1 , generic
2745	write (*,fmt='(i2,2x,f12.3,1x,g15.8)' ) &
2746	k,po(2,k,2),fo(2,k,2)
2747	end do
2748	print *,' '
2749
2750	call vert_interp ( fo , po , fn_interp , pn , &
2751	generic , 'T' , &
2752	interp_type , lagrange_order , &
2753	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
2754	zap_close_levels , force_sfc_in_vinterp , &
2755	ids , ide , jds , jde , kds , kde , &
2756	ims , ime , jms , jme , kms , kme , &
2757	its , ite , jts , jte , kts , kte )
2758
2759	print *,'Multi-Order Interpolator'
2760	print *,'------------------------------------'
2761	print *,' '
2762	print *,'Level Pressure Field Field Field'
2763	print *,' (Pa) Calc Interp Diff'
2764	print *,'------------------------------------'
2765	print *,' '
2766	do k = kts , kte-1
2767	write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) &
2768	k,pn(2,k,2),fn_calc(2,k,2),fn_interp(2,k,2),fn_calc(2,k,2)-fn_interp(2,k,2)
2769	end do
2770
2771	call vert_interp_old ( fo , po , fn_interp , pn , &
2772	generic , 'T' , &
2773	interp_type , lagrange_order , &
2774	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
2775	zap_close_levels , force_sfc_in_vinterp , &
2776	ids , ide , jds , jde , kds , kde , &
2777	ims , ime , jms , jme , kms , kme , &
2778	its , ite , jts , jte , kts , kte )
2779
2780	print *,'Linear Interpolator'
2781	print *,'------------------------------------'
2782	print *,' '
2783	print *,'Level Pressure Field Field Field'
2784	print *,' (Pa) Calc Interp Diff'
2785	print *,'------------------------------------'
2786	print *,' '
2787	do k = kts , kte-1
2788	write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) &
2789	k,pn(2,k,2),fn_calc(2,k,2),fn_interp(2,k,2),fn_calc(2,k,2)-fn_interp(2,k,2)
2790	end do
2791	end do
2792
2793	end program vint
2794
2795	subroutine wrf_error_fatal (string)
2796	character (len=*) :: string
2797	print *,string
2798	stop
2799	end subroutine wrf_error_fatal
2800
2801	subroutine fillitup ( fo , po , fn , pn , &
2802	ids , ide , jds , jde , kds , kde , &
2803	ims , ime , jms , jme , kms , kme , &
2804	its , ite , jts , jte , kts , kte , &
2805	generic , lagrange_order )
2806
2807	implicit none
2808
2809	integer , intent(in) :: ids , ide , jds , jde , kds , kde , &
2810	ims , ime , jms , jme , kms , kme , &
2811	its , ite , jts , jte , kts , kte
2812
2813	integer , intent(in) :: generic , lagrange_order
2814
2815	real , dimension(ims:ime,generic,jms:jme) , intent(out) :: fo , po
2816	real , dimension(ims:ime,kms:kme,jms:jme) , intent(out) :: fn , pn
2817
2818	integer :: i , j , k
2819
2820	real , parameter :: piov2 = 3.14159265358 / 2.
2821
2822	k = 1
2823	do j = jts , jte
2824	do i = its , ite
2825	po(i,k,j) = 102000.
2826	end do
2827	end do
2828
2829	do k = 2 , generic
2830	do j = jts , jte
2831	do i = its , ite
2832	po(i,k,j) = ( 5000. * ( 1 - (k-1) ) + 100000. * ( (k-1) - (generic-1) ) ) / (1. - real(generic-1) )
2833	end do
2834	end do
2835	end do
2836
2837	if ( lagrange_order .eq. 1 ) then
2838	do k = 1 , generic
2839	do j = jts , jte
2840	do i = its , ite
2841	fo(i,k,j) = po(i,k,j)
2842	! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. )
2843	end do
2844	end do
2845	end do
2846	else if ( lagrange_order .eq. 2 ) then
2847	do k = 1 , generic
2848	do j = jts , jte
2849	do i = its , ite
2850	fo(i,k,j) = (((po(i,k,j)-5000.)/102000.)((102000.-po(i,k,j))/102000.))102000.
2851	! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. )
2852	end do
2853	end do
2854	end do
2855	end if
2856
2857	!!!!!!!!!!!!
2858
2859	do k = kts , kte
2860	do j = jts , jte
2861	do i = its , ite
2862	pn(i,k,j) = ( 5000. * ( 0 - (k-1) ) + 102000. * ( (k-1) - (kte-1) ) ) / (-1. * real(kte-1) )
2863	end do
2864	end do
2865	end do
2866
2867	do k = kts , kte-1
2868	do j = jts , jte
2869	do i = its , ite
2870	pn(i,k,j) = ( pn(i,k,j) + pn(i,k+1,j) ) /2.
2871	end do
2872	end do
2873	end do
2874
2875
2876	if ( lagrange_order .eq. 1 ) then
2877	do k = kts , kte-1
2878	do j = jts , jte
2879	do i = its , ite
2880	fn(i,k,j) = pn(i,k,j)
2881	! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. )
2882	end do
2883	end do
2884	end do
2885	else if ( lagrange_order .eq. 2 ) then
2886	do k = kts , kte-1
2887	do j = jts , jte
2888	do i = its , ite
2889	fn(i,k,j) = (((pn(i,k,j)-5000.)/102000.)((102000.-pn(i,k,j))/102000.))102000.
2890	! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. )
2891	end do
2892	end do
2893	end do
2894	end if
2895
2896	end subroutine fillitup
2897
2898	#endif
2899
2900	!---------------------------------------------------------------------
2901
2902	SUBROUTINE vert_interp ( fo , po , fnew , pnu , &
2903	generic , var_type , &
2904	interp_type , lagrange_order , extrap_type , &
2905	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
2906	zap_close_levels , force_sfc_in_vinterp , &
2907	ids , ide , jds , jde , kds , kde , &
2908	ims , ime , jms , jme , kms , kme , &
2909	its , ite , jts , jte , kts , kte )
2910
2911	! Vertically interpolate the new field. The original field on the original
2912	! pressure levels is provided, and the new pressure surfaces to interpolate to.
2913
2914	IMPLICIT NONE
2915
2916	INTEGER , INTENT(IN) :: interp_type , lagrange_order , extrap_type
2917	LOGICAL , INTENT(IN) :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
2918	REAL , INTENT(IN) :: zap_close_levels
2919	INTEGER , INTENT(IN) :: force_sfc_in_vinterp
2920	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2921	ims , ime , jms , jme , kms , kme , &
2922	its , ite , jts , jte , kts , kte
2923	INTEGER , INTENT(IN) :: generic
2924
2925	CHARACTER (LEN=1) :: var_type
2926
2927	REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: fo , po
2928	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu
2929	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew
2930
2931	REAL , DIMENSION(ims:ime,generic,jms:jme) :: forig , porig
2932	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew
2933
2934	! Local vars
2935
2936	INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2 , knext
2937	INTEGER :: istart , iend , jstart , jend , kstart , kend
2938	INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below
2939	INTEGER , DIMENSION(ims:ime ) :: ks
2940	INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc
2941	INTEGER :: count , zap , zap_below , zap_above , kst , kcount
2942	INTEGER :: kinterp_start , kinterp_end , sfc_level
2943
2944	LOGICAL :: any_below_ground
2945
2946	REAL :: p1 , p2 , pn, hold
2947	REAL , DIMENSION(1:generic) :: ordered_porig , ordered_forig
2948	REAL , DIMENSION(kts:kte) :: ordered_pnew , ordered_fnew
2949
2950	! Excluded middle.
2951
2952	LOGICAL :: any_valid_points
2953	INTEGER :: i_valid , j_valid
2954	LOGICAL :: flip_data_required
2955
2956	! Horiontal loop bounds for different variable types.
2957
2958	IF ( var_type .EQ. 'U' ) THEN
2959	istart = its
2960	iend = ite
2961	jstart = jts
2962	jend = MIN(jde-1,jte)
2963	kstart = kts
2964	kend = kte-1
2965	DO j = jstart,jend
2966	DO k = 1,generic
2967	DO i = MAX(ids+1,its) , MIN(ide-1,ite)
2968	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2969	porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5
2970	END DO
2971	END DO
2972	IF ( ids .EQ. its ) THEN
2973	DO k = 1,generic
2974	porig(its,k,j) = po(its,k,j)
2975	END DO
2976	END IF
2977	IF ( ide .EQ. ite ) THEN
2978	DO k = 1,generic
2979	porig(ite,k,j) = po(ite-1,k,j)
2980	END DO
2981	END IF
2982
2983	DO k = kstart,kend
2984	DO i = MAX(ids+1,its) , MIN(ide-1,ite)
2985	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2986	pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5
2987	END DO
2988	END DO
2989	IF ( ids .EQ. its ) THEN
2990	DO k = kstart,kend
2991	pnew(its,k,j) = pnu(its,k,j)
2992	END DO
2993	END IF
2994	IF ( ide .EQ. ite ) THEN
2995	DO k = kstart,kend
2996	pnew(ite,k,j) = pnu(ite-1,k,j)
2997	END DO
2998	END IF
2999	END DO
3000	ELSE IF ( var_type .EQ. 'V' ) THEN
3001	istart = its
3002	iend = MIN(ide-1,ite)
3003	jstart = jts
3004	jend = jte
3005	kstart = kts
3006	kend = kte-1
3007	DO i = istart,iend
3008	DO k = 1,generic
3009	DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
3010	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3011	porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5
3012	END DO
3013	END DO
3014	IF ( jds .EQ. jts ) THEN
3015	DO k = 1,generic
3016	porig(i,k,jts) = po(i,k,jts)
3017	END DO
3018	END IF
3019	IF ( jde .EQ. jte ) THEN
3020	DO k = 1,generic
3021	porig(i,k,jte) = po(i,k,jte-1)
3022	END DO
3023	END IF
3024
3025	DO k = kstart,kend
3026	DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
3027	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3028	pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5
3029	END DO
3030	END DO
3031	IF ( jds .EQ. jts ) THEN
3032	DO k = kstart,kend
3033	pnew(i,k,jts) = pnu(i,k,jts)
3034	END DO
3035	END IF
3036	IF ( jde .EQ. jte ) THEN
3037	DO k = kstart,kend
3038	pnew(i,k,jte) = pnu(i,k,jte-1)
3039	END DO
3040	END IF
3041	END DO
3042	ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN
3043	istart = its
3044	iend = MIN(ide-1,ite)
3045	jstart = jts
3046	jend = MIN(jde-1,jte)
3047	kstart = kts
3048	kend = kte
3049	DO j = jstart,jend
3050	DO k = 1,generic
3051	DO i = istart,iend
3052	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3053	porig(i,k,j) = po(i,k,j)
3054	END DO
3055	END DO
3056
3057	DO k = kstart,kend
3058	DO i = istart,iend
3059	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3060	pnew(i,k,j) = pnu(i,k,j)
3061	END DO
3062	END DO
3063	END DO
3064	ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN
3065	istart = its
3066	iend = MIN(ide-1,ite)
3067	jstart = jts
3068	jend = MIN(jde-1,jte)
3069	kstart = kts
3070	kend = kte-1
3071	DO j = jstart,jend
3072	DO k = 1,generic
3073	DO i = istart,iend
3074	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3075	porig(i,k,j) = po(i,k,j)
3076	END DO
3077	END DO
3078
3079	DO k = kstart,kend
3080	DO i = istart,iend
3081	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3082	pnew(i,k,j) = pnu(i,k,j)
3083	END DO
3084	END DO
3085	END DO
3086	ELSE
3087	istart = its
3088	iend = MIN(ide-1,ite)
3089	jstart = jts
3090	jend = MIN(jde-1,jte)
3091	kstart = kts
3092	kend = kte-1
3093	DO j = jstart,jend
3094	DO k = 1,generic
3095	DO i = istart,iend
3096	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3097	porig(i,k,j) = po(i,k,j)
3098	END DO
3099	END DO
3100
3101	DO k = kstart,kend
3102	DO i = istart,iend
3103	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3104	pnew(i,k,j) = pnu(i,k,j)
3105	END DO
3106	END DO
3107	END DO
3108	END IF
3109
3110	! We need to find if there are any valid non-excluded-middle points in this
3111	! tile. If so, then we need to hang on to a valid i,j location.
3112
3113	any_valid_points = .false.
3114	find_valid : DO j = jstart , jend
3115	DO i = istart , iend
3116	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3117	any_valid_points = .true.
3118	i_valid = i
3119	j_valid = j
3120	EXIT find_valid
3121	END DO
3122	END DO find_valid
3123	IF ( .NOT. any_valid_points ) THEN
3124	RETURN
3125	END IF
3126
3127	IF ( porig(i_valid,2,j_valid) .LT. porig(i_valid,generic,j_valid) ) THEN
3128	flip_data_required = .true.
3129	ELSE
3130	flip_data_required = .false.
3131	END IF
3132
3133	DO j = jstart , jend
3134
3135	! The lowest level is the surface. Levels 2 through "generic" are supposed to
3136	! be "bottom-up". Flip if they are not. This is based on the input pressure
3137	! array.
3138
3139	IF ( flip_data_required ) THEN
3140	DO kn = 2 , ( generic + 1 ) / 2
3141	DO i = istart , iend
3142	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3143	hold = porig(i,kn,j)
3144	porig(i,kn,j) = porig(i,generic+2-kn,j)
3145	porig(i,generic+2-kn,j) = hold
3146	forig(i,kn,j) = fo (i,generic+2-kn,j)
3147	forig(i,generic+2-kn,j) = fo (i,kn,j)
3148	END DO
3149	END DO
3150	DO i = istart , iend
3151	forig(i,1,j) = fo (i,1,j)
3152	END DO
3153	IF ( MOD(generic,2) .EQ. 0 ) THEN
3154	k=generic/2 + 1
3155	DO i = istart , iend
3156	forig(i,k,j) = fo (i,k,j)
3157	END DO
3158	END IF
3159	ELSE
3160	DO kn = 1 , generic
3161	DO i = istart , iend
3162	forig(i,kn,j) = fo (i,kn,j)
3163	END DO
3164	END DO
3165	END IF
3166
3167	! Skip all of the levels below ground in the original data based upon the surface pressure.
3168	! The ko_above_sfc is the index in the pressure array that is above the surface. If there
3169	! are no levels underground, this is index = 2. The remaining levels are eligible for use
3170	! in the vertical interpolation.
3171
3172	DO i = istart , iend
3173	ko_above_sfc(i) = -1
3174	END DO
3175	DO ko = kstart+1 , generic
3176	DO i = istart , iend
3177	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3178	IF ( ko_above_sfc(i) .EQ. -1 ) THEN
3179	IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN
3180	ko_above_sfc(i) = ko
3181	END IF
3182	END IF
3183	END DO
3184	END DO
3185
3186	! Piece together columns of the original input data. Pass the vertical columns to
3187	! the iterpolator.
3188
3189	DO i = istart , iend
3190	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
3191
3192	! If the surface value is in the middle of the array, three steps: 1) do the
3193	! values below the ground (this is just to catch the occasional value that is
3194	! inconsistently below the surface based on input data), 2) do the surface level, then
3195	! 3) add in the levels that are above the surface. For the levels next to the surface,
3196	! we check to remove any levels that are "too close". When building the column of input
3197	! pressures, we also attend to the request for forcing the surface analysis to be used
3198	! in a few lower eta-levels.
3199
3200	! Fill in the column from up to the level just below the surface with the input
3201	! presssure and the input field (orig or old, which ever). For an isobaric input
3202	! file, this data is isobaric.
3203
3204	! How many levels have we skipped in the input column.
3205
3206	zap = 0
3207	zap_below = 0
3208	zap_above = 0
3209
3210	IF ( ko_above_sfc(i) .GT. 2 ) THEN
3211	count = 1
3212	DO ko = 2 , ko_above_sfc(i)-1
3213	ordered_porig(count) = porig(i,ko,j)
3214	ordered_forig(count) = forig(i,ko,j)
3215	count = count + 1
3216	END DO
3217
3218	! Make sure the pressure just below the surface is not "too close", this
3219	! will cause havoc with the higher order interpolators. In case of a "too close"
3220	! instance, we toss out the offending level (NOT the surface one) by simply
3221	! decrementing the accumulating loop counter.
3222
3223	IF ( ordered_porig(count-1) - porig(i,1,j) .LT. zap_close_levels ) THEN
3224	count = count -1
3225	zap = 1
3226	zap_below = 1
3227	END IF
3228
3229	! Add in the surface values.
3230
3231	ordered_porig(count) = porig(i,1,j)
3232	ordered_forig(count) = forig(i,1,j)
3233	count = count + 1
3234
3235	! A usual way to do the vertical interpolation is to pay more attention to the
3236	! surface data. Why? Well it has about 20x the density as the upper air, so we
3237	! hope the analysis is better there. We more strongly use this data by artificially
3238	! tossing out levels above the surface that are beneath a certain number of prescribed
3239	! eta levels at this (i,j). The "zap" value is how many levels of input we are
3240	! removing, which is used to tell the interpolator how many valid values are in
3241	! the column. The "count" value is the increment to the index of levels, and is
3242	! only used for assignments.
3243
3244	IF ( force_sfc_in_vinterp .GT. 0 ) THEN
3245
3246	! Get the pressure at the eta level. We want to remove all input pressure levels
3247	! between the level above the surface to the pressure at this eta surface. That
3248	! forces the surface value to be used through the selected eta level. Keep track
3249	! of two things: the level to use above the eta levels, and how many levels we are
3250	! skipping.
3251
3252	knext = ko_above_sfc(i)
3253	find_level : DO ko = ko_above_sfc(i) , generic
3254	IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN
3255	knext = ko
3256	exit find_level
3257	ELSE
3258	zap = zap + 1
3259	zap_above = zap_above + 1
3260	END IF
3261	END DO find_level
3262
3263	! No request for special interpolation, so we just assign the next level to use
3264	! above the surface as, ta da, the first level above the surface. I know, wow.
3265
3266	ELSE
3267	knext = ko_above_sfc(i)
3268	END IF
3269
3270	! One more time, make sure the pressure just above the surface is not "too close", this
3271	! will cause havoc with the higher order interpolators. In case of a "too close"
3272	! instance, we toss out the offending level above the surface (NOT the surface one) by simply
3273	! incrementing the loop counter. Here, count-1 is the surface level and knext is either
3274	! the next level up OR it is the level above the prescribed number of eta surfaces.
3275
3276	IF ( ordered_porig(count-1) - porig(i,knext,j) .LT. zap_close_levels ) THEN
3277	kst = knext+1
3278	zap = zap + 1
3279	zap_above = zap_above + 1
3280	ELSE
3281	kst = knext
3282	END IF
3283
3284	DO ko = kst , generic
3285	ordered_porig(count) = porig(i,ko,j)
3286	ordered_forig(count) = forig(i,ko,j)
3287	count = count + 1
3288	END DO
3289
3290	! This is easy, the surface is the lowest level, just stick them in, in this order. OK,
3291	! there are a couple of subtleties. We have to check for that special interpolation that
3292	! skips some input levels so that the surface is used for the lowest few eta levels. Also,
3293	! we must make sure that we still do not have levels that are "too close" together.
3294
3295	ELSE
3296
3297	! Initialize no input levels have yet been removed from consideration.
3298
3299	zap = 0
3300
3301	! The surface is the lowest level, so it gets set right away to location 1.
3302
3303	ordered_porig(1) = porig(i,1,j)
3304	ordered_forig(1) = forig(i,1,j)
3305
3306	! We start filling in the array at loc 2, as in just above the level we just stored.
3307
3308	count = 2
3309
3310	! Are we forcing the interpolator to skip valid input levels so that the
3311	! surface data is used through more levels? Essentially as above.
3312
3313	IF ( force_sfc_in_vinterp .GT. 0 ) THEN
3314	knext = 2
3315	find_level2: DO ko = 2 , generic
3316	IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN
3317	knext = ko
3318	exit find_level2
3319	ELSE
3320	zap = zap + 1
3321	zap_above = zap_above + 1
3322	END IF
3323	END DO find_level2
3324	ELSE
3325	knext = 2
3326	END IF
3327
3328	! Fill in the data above the surface. The "knext" index is either the one
3329	! just above the surface OR it is the index associated with the level that
3330	! is just above the pressure at this (i,j) of the top eta level that is to
3331	! be directly impacted with the surface level in interpolation.
3332
3333	DO ko = knext , generic
3334	IF ( ( ordered_porig(count-1) - porig(i,ko,j) .LT. zap_close_levels ) .AND. &
3335	( ko .LT. generic ) ) THEN
3336	zap = zap + 1
3337	zap_above = zap_above + 1
3338	CYCLE
3339	END IF
3340	ordered_porig(count) = porig(i,ko,j)
3341	ordered_forig(count) = forig(i,ko,j)
3342	count = count + 1
3343	END DO
3344
3345	END IF
3346
3347	! Now get the column of the "new" pressure data. So, this one is easy.
3348
3349	DO kn = kstart , kend
3350	ordered_pnew(kn) = pnew(i,kn,j)
3351	END DO
3352
3353	! How many levels (count) are we shipping to the Lagrange interpolator.
3354
3355	IF ( ( use_levels_below_ground ) .AND. ( use_surface ) ) THEN
3356
3357	! Use all levels, including the input surface, and including the pressure
3358	! levels below ground. We know to stop when we have reached the top of
3359	! the input pressure data.
3360
3361	count = 0
3362	find_how_many_1 : DO ko = 1 , generic
3363	IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN
3364	count = count + 1
3365	EXIT find_how_many_1
3366	ELSE
3367	count = count + 1
3368	END IF
3369	END DO find_how_many_1
3370	kinterp_start = 1
3371	kinterp_end = kinterp_start + count - 1
3372
3373	ELSE IF ( ( use_levels_below_ground ) .AND. ( .NOT. use_surface ) ) THEN
3374
3375	! Use all levels (excluding the input surface) and including the pressure
3376	! levels below ground. We know to stop when we have reached the top of
3377	! the input pressure data.
3378
3379	count = 0
3380	find_sfc_2 : DO ko = 1 , generic
3381	IF ( porig(i,1,j) .EQ. ordered_porig(ko) ) THEN
3382	sfc_level = ko
3383	EXIT find_sfc_2
3384	END IF
3385	END DO find_sfc_2
3386
3387	DO ko = sfc_level , generic-1
3388	ordered_porig(ko) = ordered_porig(ko+1)
3389	ordered_forig(ko) = ordered_forig(ko+1)
3390	END DO
3391	ordered_porig(generic) = 1.E-5
3392	ordered_forig(generic) = 1.E10
3393
3394	count = 0
3395	find_how_many_2 : DO ko = 1 , generic
3396	IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN
3397	count = count + 1
3398	EXIT find_how_many_2
3399	ELSE
3400	count = count + 1
3401	END IF
3402	END DO find_how_many_2
3403	kinterp_start = 1
3404	kinterp_end = kinterp_start + count - 1
3405
3406	ELSE IF ( ( .NOT. use_levels_below_ground ) .AND. ( use_surface ) ) THEN
3407
3408	! Use all levels above the input surface pressure.
3409
3410	kcount = ko_above_sfc(i)-1-zap_below
3411	count = 0
3412	DO ko = 1 , generic
3413	IF ( porig(i,ko,j) .EQ. ordered_porig(kcount) ) THEN
3414	! write (6,fmt='(f11.3,f11.3,g11.5)') porig(i,ko,j),ordered_porig(kcount),ordered_forig(kcount)
3415	kcount = kcount + 1
3416	count = count + 1
3417	ELSE
3418	! write (6,fmt='(f11.3 )') porig(i,ko,j)
3419	END IF
3420	END DO
3421	kinterp_start = ko_above_sfc(i)-1-zap_below
3422	kinterp_end = kinterp_start + count - 1
3423
3424	END IF
3425
3426	! The polynomials are either in pressure or LOG(pressure).
3427
3428	IF ( interp_type .EQ. 1 ) THEN
3429	CALL lagrange_setup ( var_type , interp_type , &
3430	ordered_porig(kinterp_start:kinterp_end) , &
3431	ordered_forig(kinterp_start:kinterp_end) , &
3432	count , lagrange_order , extrap_type , &
3433	ordered_pnew(kstart:kend) , ordered_fnew , kend-kstart+1 ,i,j)
3434	ELSE
3435	CALL lagrange_setup ( var_type , interp_type , &
3436	LOG(ordered_porig(kinterp_start:kinterp_end)) , &
3437	ordered_forig(kinterp_start:kinterp_end) , &
3438	count , lagrange_order , extrap_type , &
3439	LOG(ordered_pnew(kstart:kend)) , ordered_fnew , kend-kstart+1 ,i,j)
3440	END IF
3441
3442	! Save the computed data.
3443
3444	DO kn = kstart , kend
3445	fnew(i,kn,j) = ordered_fnew(kn)
3446	END DO
3447
3448	! There may have been a request to have the surface data from the input field
3449	! to be assigned as to the lowest eta level. This assumes thin layers (usually
3450	! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V).
3451
3452	IF ( lowest_lev_from_sfc ) THEN
3453	fnew(i,1,j) = forig(i,ko_above_sfc(i)-1,j)
3454	END IF
3455
3456	END DO
3457
3458	END DO
3459
3460	END SUBROUTINE vert_interp
3461
3462	!---------------------------------------------------------------------
3463
3464	SUBROUTINE vert_interp_old ( forig , po , fnew , pnu , &
3465	generic , var_type , &
3466	interp_type , lagrange_order , extrap_type , &
3467	lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
3468	zap_close_levels , force_sfc_in_vinterp , &
3469	ids , ide , jds , jde , kds , kde , &
3470	ims , ime , jms , jme , kms , kme , &
3471	its , ite , jts , jte , kts , kte )
3472
3473	! Vertically interpolate the new field. The original field on the original
3474	! pressure levels is provided, and the new pressure surfaces to interpolate to.
3475
3476	IMPLICIT NONE
3477
3478	INTEGER , INTENT(IN) :: interp_type , lagrange_order , extrap_type
3479	LOGICAL , INTENT(IN) :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
3480	REAL , INTENT(IN) :: zap_close_levels
3481	INTEGER , INTENT(IN) :: force_sfc_in_vinterp
3482	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3483	ims , ime , jms , jme , kms , kme , &
3484	its , ite , jts , jte , kts , kte
3485	INTEGER , INTENT(IN) :: generic
3486
3487	CHARACTER (LEN=1) :: var_type
3488
3489	REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: forig , po
3490	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu
3491	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew
3492
3493	REAL , DIMENSION(ims:ime,generic,jms:jme) :: porig
3494	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew
3495
3496	! Local vars
3497
3498	INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2
3499	INTEGER :: istart , iend , jstart , jend , kstart , kend
3500	INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below
3501	INTEGER , DIMENSION(ims:ime ) :: ks
3502	INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc
3503
3504	LOGICAL :: any_below_ground
3505
3506	REAL :: p1 , p2 , pn
3507	integer vert_extrap
3508	vert_extrap = 0
3509
3510	! Horiontal loop bounds for different variable types.
3511
3512	IF ( var_type .EQ. 'U' ) THEN
3513	istart = its
3514	iend = ite
3515	jstart = jts
3516	jend = MIN(jde-1,jte)
3517	kstart = kts
3518	kend = kte-1
3519	DO j = jstart,jend
3520	DO k = 1,generic
3521	DO i = MAX(ids+1,its) , MIN(ide-1,ite)
3522	porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5
3523	END DO
3524	END DO
3525	IF ( ids .EQ. its ) THEN
3526	DO k = 1,generic
3527	porig(its,k,j) = po(its,k,j)
3528	END DO
3529	END IF
3530	IF ( ide .EQ. ite ) THEN
3531	DO k = 1,generic
3532	porig(ite,k,j) = po(ite-1,k,j)
3533	END DO
3534	END IF
3535
3536	DO k = kstart,kend
3537	DO i = MAX(ids+1,its) , MIN(ide-1,ite)
3538	pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5
3539	END DO
3540	END DO
3541	IF ( ids .EQ. its ) THEN
3542	DO k = kstart,kend
3543	pnew(its,k,j) = pnu(its,k,j)
3544	END DO
3545	END IF
3546	IF ( ide .EQ. ite ) THEN
3547	DO k = kstart,kend
3548	pnew(ite,k,j) = pnu(ite-1,k,j)
3549	END DO
3550	END IF
3551	END DO
3552	ELSE IF ( var_type .EQ. 'V' ) THEN
3553	istart = its
3554	iend = MIN(ide-1,ite)
3555	jstart = jts
3556	jend = jte
3557	kstart = kts
3558	kend = kte-1
3559	DO i = istart,iend
3560	DO k = 1,generic
3561	DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
3562	porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5
3563	END DO
3564	END DO
3565	IF ( jds .EQ. jts ) THEN
3566	DO k = 1,generic
3567	porig(i,k,jts) = po(i,k,jts)
3568	END DO
3569	END IF
3570	IF ( jde .EQ. jte ) THEN
3571	DO k = 1,generic
3572	porig(i,k,jte) = po(i,k,jte-1)
3573	END DO
3574	END IF
3575
3576	DO k = kstart,kend
3577	DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
3578	pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5
3579	END DO
3580	END DO
3581	IF ( jds .EQ. jts ) THEN
3582	DO k = kstart,kend
3583	pnew(i,k,jts) = pnu(i,k,jts)
3584	END DO
3585	END IF
3586	IF ( jde .EQ. jte ) THEN
3587	DO k = kstart,kend
3588	pnew(i,k,jte) = pnu(i,k,jte-1)
3589	END DO
3590	END IF
3591	END DO
3592	ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN
3593	istart = its
3594	iend = MIN(ide-1,ite)
3595	jstart = jts
3596	jend = MIN(jde-1,jte)
3597	kstart = kts
3598	kend = kte
3599	DO j = jstart,jend
3600	DO k = 1,generic
3601	DO i = istart,iend
3602	porig(i,k,j) = po(i,k,j)
3603	END DO
3604	END DO
3605
3606	DO k = kstart,kend
3607	DO i = istart,iend
3608	pnew(i,k,j) = pnu(i,k,j)
3609	END DO
3610	END DO
3611	END DO
3612	ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN
3613	istart = its
3614	iend = MIN(ide-1,ite)
3615	jstart = jts
3616	jend = MIN(jde-1,jte)
3617	kstart = kts
3618	kend = kte-1
3619	DO j = jstart,jend
3620	DO k = 1,generic
3621	DO i = istart,iend
3622	porig(i,k,j) = po(i,k,j)
3623	END DO
3624	END DO
3625
3626	DO k = kstart,kend
3627	DO i = istart,iend
3628	pnew(i,k,j) = pnu(i,k,j)
3629	END DO
3630	END DO
3631	END DO
3632	ELSE
3633	istart = its
3634	iend = MIN(ide-1,ite)
3635	jstart = jts
3636	jend = MIN(jde-1,jte)
3637	kstart = kts
3638	kend = kte-1
3639	DO j = jstart,jend
3640	DO k = 1,generic
3641	DO i = istart,iend
3642	porig(i,k,j) = po(i,k,j)
3643	END DO
3644	END DO
3645
3646	DO k = kstart,kend
3647	DO i = istart,iend
3648	pnew(i,k,j) = pnu(i,k,j)
3649	END DO
3650	END DO
3651	END DO
3652	END IF
3653
3654	DO j = jstart , jend
3655
3656	! Skip all of the levels below ground in the original data based upon the surface pressure.
3657	! The ko_above_sfc is the index in the pressure array that is above the surface. If there
3658	! are no levels underground, this is index = 2. The remaining levels are eligible for use
3659	! in the vertical interpolation.
3660
3661	DO i = istart , iend
3662	ko_above_sfc(i) = -1
3663	END DO
3664	DO ko = kstart+1 , kend
3665	DO i = istart , iend
3666	IF ( ko_above_sfc(i) .EQ. -1 ) THEN
3667	IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN
3668	ko_above_sfc(i) = ko
3669	END IF
3670	END IF
3671	END DO
3672	END DO
3673
3674	! Initialize interpolation location. These are the levels in the original pressure
3675	! data that are physically below and above the targeted new pressure level.
3676
3677	DO kn = kts , kte
3678	DO i = its , ite
3679	k_above(i,kn) = -1
3680	k_below(i,kn) = -2
3681	END DO
3682	END DO
3683
3684	! Starting location is no lower than previous found location. This is for O(n logn)
3685	! and not O(n^2), where n is the number of vertical levels to search.
3686
3687	DO i = its , ite
3688	ks(i) = 1
3689	END DO
3690
3691	! Find trapping layer for interpolation. The kn index runs through all of the "new"
3692	! levels of data.
3693
3694	DO kn = kstart , kend
3695
3696	DO i = istart , iend
3697
3698	! For each "new" level (kn), we search to find the trapping levels in the "orig"
3699	! data. Most of the time, the "new" levels are the eta surfaces, and the "orig"
3700	! levels are the input pressure levels.
3701
3702	found_trap_above : DO ko = ks(i) , generic-1
3703
3704	! Because we can have levels in the interpolation that are not valid,
3705	! let's toss out any candidate orig pressure values that are below ground
3706	! based on the surface pressure. If the level =1, then this IS the surface
3707	! level, so we HAVE to keep that one, but maybe not the ones above. If the
3708	! level (ks) is NOT=1, then we have to just CYCLE our loop to find a legit
3709	! below-pressure value. If we are not below ground, then we choose two
3710	! neighboring levels to test whether they surround the new pressure level.
3711
3712	! The input trapping levels that we are trying is the surface and the first valid
3713	! level above the surface.
3714
3715	IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .EQ. 1 ) ) THEN
3716	ko_1 = ko
3717	ko_2 = ko_above_sfc(i)
3718
3719	! The "below" level is underground, cycle until we get to a valid pressure
3720	! above ground.
3721
3722	ELSE IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .NE. 1 ) ) THEN
3723	CYCLE found_trap_above
3724
3725	! The "below" level is above the surface, so we are in the clear to test these
3726	! two levels out.
3727
3728	ELSE
3729	ko_1 = ko
3730	ko_2 = ko+1
3731
3732	END IF
3733
3734	! The test of the candidate levels: "below" has to have a larger pressure, and
3735	! "above" has to have a smaller pressure.
3736
3737	! OK, we found the correct two surrounding levels. The locations are saved for use in the
3738	! interpolation.
3739
3740	IF ( ( porig(i,ko_1,j) .GE. pnew(i,kn,j) ) .AND. &
3741	( porig(i,ko_2,j) .LT. pnew(i,kn,j) ) ) THEN
3742	k_above(i,kn) = ko_2
3743	k_below(i,kn) = ko_1
3744	ks(i) = ko_1
3745	EXIT found_trap_above
3746
3747	! What do we do is we need to extrapolate the data underground? This happens when the
3748	! lowest pressure that we have is physically "above" the new target pressure. Our
3749	! actions depend on the type of variable we are interpolating.
3750
3751	ELSE IF ( porig(i,1,j) .LT. pnew(i,kn,j) ) THEN
3752
3753	! For horizontal winds and moisture, we keep a constant value under ground.
3754
3755	IF ( ( var_type .EQ. 'U' ) .OR. &
3756	( var_type .EQ. 'V' ) .OR. &
3757	( var_type .EQ. 'Q' ) ) THEN
3758	k_above(i,kn) = 1
3759	ks(i) = 1
3760
3761	! For temperature and height, we extrapolate the data. Hopefully, we are not
3762	! extrapolating too far. For pressure level input, the eta levels are always
3763	! contained within the surface to p_top levels, so no extrapolation is ever
3764	! required.
3765
3766	ELSE IF ( ( var_type .EQ. 'Z' ) .OR. &
3767	( var_type .EQ. 'T' ) ) THEN
3768	k_above(i,kn) = ko_above_sfc(i)
3769	k_below(i,kn) = 1
3770	ks(i) = 1
3771
3772	! Just a catch all right now.
3773
3774	ELSE
3775	k_above(i,kn) = 1
3776	ks(i) = 1
3777	END IF
3778
3779	EXIT found_trap_above
3780
3781	! The other extrapolation that might be required is when we are going above the
3782	! top level of the input data. Usually this means we chose a P_PTOP value that
3783	! was inappropriate, and we should stop and let someone fix this mess.
3784
3785	ELSE IF ( porig(i,generic,j) .GT. pnew(i,kn,j) ) THEN
3786	print *,'data is too high, try a lower p_top'
3787	print *,'pnew=',pnew(i,kn,j)
3788	print *,'porig=',porig(i,:,j)
3789	CALL wrf_error_fatal ('requested p_top is higher than input data, lower p_top')
3790
3791	END IF
3792	END DO found_trap_above
3793	END DO
3794	END DO
3795
3796	! Linear vertical interpolation.
3797
3798	DO kn = kstart , kend
3799	DO i = istart , iend
3800	IF ( k_above(i,kn) .EQ. 1 ) THEN
3801	fnew(i,kn,j) = forig(i,1,j)
3802	ELSE
3803	k2 = MAX ( k_above(i,kn) , 2)
3804	k1 = MAX ( k_below(i,kn) , 1)
3805	IF ( k1 .EQ. k2 ) THEN
3806	CALL wrf_error_fatal ( 'identical values in the interp, bad for divisions' )
3807	END IF
3808	IF ( interp_type .EQ. 1 ) THEN
3809	p1 = porig(i,k1,j)
3810	p2 = porig(i,k2,j)
3811	pn = pnew(i,kn,j)
3812	ELSE IF ( interp_type .EQ. 2 ) THEN
3813	p1 = ALOG(porig(i,k1,j))
3814	p2 = ALOG(porig(i,k2,j))
3815	pn = ALOG(pnew(i,kn,j))
3816	END IF
3817	IF ( ( p1-pn) * (p2-pn) > 0. ) THEN
3818	! CALL wrf_error_fatal ( 'both trapping pressures are on the same side of the new pressure' )
3819	! CALL wrf_debug ( 0 , 'both trapping pressures are on the same side of the new pressure' )
3820	vert_extrap = vert_extrap + 1
3821	END IF
3822	fnew(i,kn,j) = ( forig(i,k1,j) * ( p2 - pn ) + &
3823	forig(i,k2,j) * ( pn - p1 ) ) / &
3824	( p2 - p1 )
3825	END IF
3826	END DO
3827	END DO
3828
3829	search_below_ground : DO kn = kstart , kend
3830	any_below_ground = .FALSE.
3831	DO i = istart , iend
3832	IF ( k_above(i,kn) .EQ. 1 ) THEN
3833	fnew(i,kn,j) = forig(i,1,j)
3834	any_below_ground = .TRUE.
3835	END IF
3836	END DO
3837	IF ( .NOT. any_below_ground ) THEN
3838	EXIT search_below_ground
3839	END IF
3840	END DO search_below_ground
3841
3842	! There may have been a request to have the surface data from the input field
3843	! to be assigned as to the lowest eta level. This assumes thin layers (usually
3844	! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V).
3845
3846	DO i = istart , iend
3847	IF ( lowest_lev_from_sfc ) THEN
3848	fnew(i,1,j) = forig(i,ko_above_sfc(i),j)
3849	END IF
3850	END DO
3851
3852	END DO
3853	print *,'VERT EXTRAP = ', vert_extrap
3854
3855	END SUBROUTINE vert_interp_old
3856
3857	!---------------------------------------------------------------------
3858
3859	SUBROUTINE lagrange_setup ( var_type , interp_type , all_x , all_y , all_dim , n , extrap_type , &
3860	target_x , target_y , target_dim ,i,j)
3861
3862	! We call a Lagrange polynomial interpolator. The parallel concerns are put off as this
3863	! is initially set up for vertical use. The purpose is an input column of pressure (all_x),
3864	! and the associated pressure level data (all_y). These are assumed to be sorted (ascending
3865	! or descending, no matter). The locations to be interpolated to are the pressures in
3866	! target_x, probably the new vertical coordinate values. The field that is output is the
3867	! target_y, which is defined at the target_x location. Mostly we expect to be 2nd order
3868	! overlapping polynomials, with only a single 2nd order method near the top and bottom.
3869	! When n=1, this is linear; when n=2, this is a second order interpolator.
3870
3871	IMPLICIT NONE
3872
3873	CHARACTER (LEN=1) :: var_type
3874	INTEGER , INTENT(IN) :: interp_type , all_dim , n , extrap_type , target_dim
3875	REAL, DIMENSION(all_dim) , INTENT(IN) :: all_x , all_y
3876	REAL , DIMENSION(target_dim) , INTENT(IN) :: target_x
3877	REAL , DIMENSION(target_dim) , INTENT(OUT) :: target_y
3878
3879	! Brought in for debug purposes, all of the computations are in a single column.
3880
3881	INTEGER , INTENT(IN) :: i,j
3882
3883	! Local vars
3884
3885	REAL , DIMENSION(n+1) :: x , y
3886	REAL :: a , b
3887	REAL :: target_y_1 , target_y_2
3888	LOGICAL :: found_loc
3889	INTEGER :: loop , loc_center_left , loc_center_right , ist , iend , target_loop
3890	INTEGER :: vboundb , vboundt
3891
3892	! Local vars for the problem of extrapolating theta below ground.
3893
3894	REAL :: temp_1 , temp_2 , temp_3 , temp_y
3895	REAL :: depth_of_extrap_in_p , avg_of_extrap_p , temp_extrap_starting_point , dhdp , dh , dt
3896	REAL , PARAMETER :: RovCp = rcp
3897	REAL , PARAMETER :: CRC_const1 = 11880.516 ! m
3898	REAL , PARAMETER :: CRC_const2 = 0.1902632 !
3899	REAL , PARAMETER :: CRC_const3 = 0.0065 ! K/km
3900	REAL, DIMENSION(all_dim) :: all_x_full
3901	REAL , DIMENSION(target_dim) :: target_x_full
3902
3903	IF ( all_dim .LT. n+1 ) THEN
3904	print *,'all_dim = ',all_dim
3905	print *,'order = ',n
3906	print *,'i,j = ',i,j
3907	print *,'p array = ',all_x
3908	print *,'f array = ',all_y
3909	print *,'p target= ',target_x
3910	CALL wrf_error_fatal ( 'troubles, the interpolating order is too large for this few input values' )
3911	END IF
3912
3913	IF ( n .LT. 1 ) THEN
3914	CALL wrf_error_fatal ( 'pal, linear is about as low as we go' )
3915	END IF
3916
3917	! We can pinch in the area of the higher order interpolation with vbound. If
3918	! vbound = 0, no pinching. If vbound = m, then we make the lower "m" and upper
3919	! "m" eta levels use a linear interpolation.
3920
3921	vboundb = 4
3922	vboundt = 0
3923
3924	! Loop over the list of target x and y values.
3925
3926	DO target_loop = 1 , target_dim
3927
3928	! Find the two trapping x values, and keep the indices.
3929
3930	found_loc = .FALSE.
3931	find_trap : DO loop = 1 , all_dim -1
3932	a = target_x(target_loop) - all_x(loop)
3933	b = target_x(target_loop) - all_x(loop+1)
3934	IF ( a*b .LE. 0.0 ) THEN
3935	loc_center_left = loop
3936	loc_center_right = loop+1
3937	found_loc = .TRUE.
3938	EXIT find_trap
3939	END IF
3940	END DO find_trap
3941
3942	IF ( ( .NOT. found_loc ) .AND. ( target_x(target_loop) .GT. all_x(1) ) ) THEN
3943
3944	! Get full pressure back so that our extrpolations make sense.
3945
3946	IF ( interp_type .EQ. 1 ) THEN
3947	all_x_full = all_x
3948	target_x_full = target_x
3949	ELSE
3950	all_x_full = EXP ( all_x )
3951	target_x_full = EXP ( target_x )
3952	END IF
3953	! Isothermal extrapolation.
3954
3955	IF ( ( extrap_type .EQ. 1 ) .AND. ( var_type .EQ. 'T' ) ) THEN
3956
3957	temp_1 = all_y(1) * ( all_x_full(1) / 100000. ) ** RovCp
3958	target_y(target_loop) = temp_1 * ( 100000. / target_x_full(target_loop) ) ** RovCp
3959
3960	! Standard atmosphere -6.5 K/km lapse rate for the extrapolation.
3961
3962	ELSE IF ( ( extrap_type .EQ. 2 ) .AND. ( var_type .EQ. 'T' ) ) THEN
3963
3964	depth_of_extrap_in_p = target_x_full(target_loop) - all_x_full(1)
3965	avg_of_extrap_p = ( target_x_full(target_loop) + all_x_full(1) ) * 0.5
3966	temp_extrap_starting_point = all_y(1) * ( all_x_full(1) / 100000. ) ** RovCp
3967	dhdp = CRC_const1 * CRC_const2 * ( avg_of_extrap_p / 100. ) ** ( CRC_const2 - 1. )
3968	dh = dhdp * ( depth_of_extrap_in_p / 100. )
3969	dt = dh * CRC_const3
3970	target_y(target_loop) = ( temp_extrap_starting_point + dt ) * ( 100000. / target_x_full(target_loop) ) ** RovCp
3971
3972	! Adiabatic extrapolation for theta.
3973
3974	ELSE IF ( ( extrap_type .EQ. 3 ) .AND. ( var_type .EQ. 'T' ) ) THEN
3975
3976	target_y(target_loop) = all_y(1)
3977
3978
3979	! Wild extrapolation for non-temperature vars.
3980
3981	ELSE IF ( extrap_type .EQ. 1 ) THEN
3982
3983	target_y(target_loop) = ( all_y(2) * ( target_x(target_loop) - all_x(3) ) + &
3984	all_y(3) * ( all_x(2) - target_x(target_loop) ) ) / &
3985	( all_x(2) - all_x(3) )
3986
3987	! Use a constant value below ground.
3988
3989	ELSE IF ( extrap_type .EQ. 2 ) THEN
3990
3991	target_y(target_loop) = all_y(1)
3992
3993	ELSE IF ( extrap_type .EQ. 3 ) THEN
3994	CALL wrf_error_fatal ( 'You are not allowed to use extrap_option #3 for any var except for theta.' )
3995
3996	END IF
3997	CYCLE
3998	ELSE IF ( .NOT. found_loc ) THEN
3999	print *,'i,j = ',i,j
4000	print *,'target pressure and value = ',target_x(target_loop),target_y(target_loop)
4001	DO loop = 1 , all_dim
4002	print *,'column of pressure and value = ',all_x(loop),all_y(loop)
4003	END DO
4004	CALL wrf_error_fatal ( 'troubles, could not find trapping x locations' )
4005	END IF
4006
4007	! Even or odd order? We can put the value in the middle if this is
4008	! an odd order interpolator. For the even guys, we'll do it twice
4009	! and shift the range one index, then get an average.
4010
4011	IF ( MOD(n,2) .NE. 0 ) THEN
4012	IF ( ( loc_center_left -(((n+1)/2)-1) .GE. 1 ) .AND. &
4013	( loc_center_right+(((n+1)/2)-1) .LE. all_dim ) ) THEN
4014	ist = loc_center_left -(((n+1)/2)-1)
4015	iend = ist + n
4016	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
4017	ELSE
4018	IF ( .NOT. found_loc ) THEN
4019	CALL wrf_error_fatal ( 'I doubt this will happen, I will only do 2nd order for now' )
4020	END IF
4021	END IF
4022
4023	ELSE IF ( ( MOD(n,2) .EQ. 0 ) .AND. &
4024	( ( target_loop .GE. 1 + vboundb ) .AND. ( target_loop .LE. target_dim - vboundt ) ) ) THEN
4025	IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. &
4026	( loc_center_right+(((n )/2) ) .LE. all_dim ) .AND. &
4027	( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. &
4028	( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN
4029	ist = loc_center_left -(((n )/2)-1)
4030	iend = ist + n
4031	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_1 )
4032	ist = loc_center_left -(((n )/2) )
4033	iend = ist + n
4034	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_2 )
4035	target_y(target_loop) = ( target_y_1 + target_y_2 ) * 0.5
4036
4037	ELSE IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. &
4038	( loc_center_right+(((n )/2) ) .LE. all_dim ) ) THEN
4039	ist = loc_center_left -(((n )/2)-1)
4040	iend = ist + n
4041	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
4042	ELSE IF ( ( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. &
4043	( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN
4044	ist = loc_center_left -(((n )/2) )
4045	iend = ist + n
4046	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
4047	ELSE
4048	CALL wrf_error_fatal ( 'unauthorized area, you should not be here' )
4049	END IF
4050
4051	ELSE IF ( MOD(n,2) .EQ. 0 ) THEN
4052	ist = loc_center_left
4053	iend = loc_center_right
4054	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , 1 , target_x(target_loop) , target_y(target_loop) )
4055
4056	END IF
4057
4058	END DO
4059
4060	END SUBROUTINE lagrange_setup
4061
4062	!---------------------------------------------------------------------
4063
4064	SUBROUTINE lagrange_interp ( x , y , n , target_x , target_y )
4065
4066	! Interpolation using Lagrange polynomials.
4067	! P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
4068	! where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
4069	! ---------------------------------------------
4070	! (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
4071
4072	IMPLICIT NONE
4073
4074	INTEGER , INTENT(IN) :: n
4075	REAL , DIMENSION(0:n) , INTENT(IN) :: x , y
4076	REAL , INTENT(IN) :: target_x
4077
4078	REAL , INTENT(OUT) :: target_y
4079
4080	! Local vars
4081
4082	INTEGER :: i , k
4083	REAL :: numer , denom , Px
4084	REAL , DIMENSION(0:n) :: Ln
4085
4086	Px = 0.
4087	DO i = 0 , n
4088	numer = 1.
4089	denom = 1.
4090	DO k = 0 , n
4091	IF ( k .EQ. i ) CYCLE
4092	numer = numer * ( target_x - x(k) )
4093	denom = denom * ( x(i) - x(k) )
4094	END DO
4095	IF ( denom .NE. 0. ) THEN
4096	Ln(i) = y(i) * numer / denom
4097	Px = Px + Ln(i)
4098	ENDIF
4099	END DO
4100	target_y = Px
4101
4102	END SUBROUTINE lagrange_interp
4103
4104	#ifndef VERT_UNIT
4105	!---------------------------------------------------------------------
4106
4107	SUBROUTINE p_dry ( mu0 , eta , pdht , pdry , full_levs , &
4108	ids , ide , jds , jde , kds , kde , &
4109	ims , ime , jms , jme , kms , kme , &
4110	its , ite , jts , jte , kts , kte )
4111
4112	! Compute reference pressure and the reference mu.
4113
4114	IMPLICIT NONE
4115
4116	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4117	ims , ime , jms , jme , kms , kme , &
4118	its , ite , jts , jte , kts , kte
4119
4120	LOGICAL :: full_levs
4121
4122	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(IN) :: mu0
4123	REAL , DIMENSION( kms:kme ) , INTENT(IN) :: eta
4124	REAL :: pdht
4125	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pdry
4126
4127	! Local vars
4128
4129	INTEGER :: i , j , k
4130	REAL , DIMENSION( kms:kme ) :: eta_h
4131
4132	IF ( full_levs ) THEN
4133	DO j = jts , MIN ( jde-1 , jte )
4134	DO k = kts , kte
4135	DO i = its , MIN (ide-1 , ite )
4136	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4137	pdry(i,k,j) = eta(k) * mu0(i,j) + pdht
4138	END DO
4139	END DO
4140	END DO
4141
4142	ELSE
4143	DO k = kts , kte-1
4144	eta_h(k) = ( eta(k) + eta(k+1) ) * 0.5
4145	END DO
4146
4147	DO j = jts , MIN ( jde-1 , jte )
4148	DO k = kts , kte-1
4149	DO i = its , MIN (ide-1 , ite )
4150	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4151	pdry(i,k,j) = eta_h(k) * mu0(i,j) + pdht
4152	END DO
4153	END DO
4154	END DO
4155	END IF
4156
4157	END SUBROUTINE p_dry
4158
4159	!---------------------------------------------------------------------
4160
4161	SUBROUTINE p_dts ( pdts , intq , psfc , p_top , &
4162	ids , ide , jds , jde , kds , kde , &
4163	ims , ime , jms , jme , kms , kme , &
4164	its , ite , jts , jte , kts , kte )
4165
4166	! Compute difference between the dry, total surface pressure and the top pressure.
4167
4168	IMPLICIT NONE
4169
4170	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4171	ims , ime , jms , jme , kms , kme , &
4172	its , ite , jts , jte , kts , kte
4173
4174	REAL , INTENT(IN) :: p_top
4175	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: psfc
4176	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: intq
4177	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(OUT) :: pdts
4178
4179	! Local vars
4180
4181	INTEGER :: i , j , k
4182
4183	DO j = jts , MIN ( jde-1 , jte )
4184	DO i = its , MIN (ide-1 , ite )
4185	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4186	pdts(i,j) = psfc(i,j) - intq(i,j) - p_top
4187	END DO
4188	END DO
4189
4190	END SUBROUTINE p_dts
4191
4192	!---------------------------------------------------------------------
4193
4194	SUBROUTINE p_dhs ( pdhs , ht , p0 , t0 , a , &
4195	ids , ide , jds , jde , kds , kde , &
4196	ims , ime , jms , jme , kms , kme , &
4197	its , ite , jts , jte , kts , kte )
4198
4199	! Compute dry, hydrostatic surface pressure.
4200
4201	IMPLICIT NONE
4202
4203	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4204	ims , ime , jms , jme , kms , kme , &
4205	its , ite , jts , jte , kts , kte
4206
4207	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(IN) :: ht
4208	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: pdhs
4209
4210	REAL , INTENT(IN) :: p0 , t0 , a
4211
4212	! Local vars
4213
4214	INTEGER :: i , j , k
4215
4216	REAL , PARAMETER :: Rd = r_d
4217
4218	DO j = jts , MIN ( jde-1 , jte )
4219	DO i = its , MIN (ide-1 , ite )
4220	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4221	pdhs(i,j) = p0 * EXP ( -t0/a + SQRT ( (t0/a)*2 - 2. g * ht(i,j)/(a * Rd) ) )
4222	END DO
4223	END DO
4224
4225	END SUBROUTINE p_dhs
4226
4227	!---------------------------------------------------------------------
4228
4229	SUBROUTINE find_p_top ( p , p_top , &
4230	ids , ide , jds , jde , kds , kde , &
4231	ims , ime , jms , jme , kms , kme , &
4232	its , ite , jts , jte , kts , kte )
4233
4234	! Find the largest pressure in the top level. This is our p_top. We are
4235	! assuming that the top level is the location where the pressure is a minimum
4236	! for each column. In cases where the top surface is not isobaric, a
4237	! communicated value must be shared in the calling routine. Also in cases
4238	! where the top surface is not isobaric, care must be taken that the new
4239	! maximum pressure is not greater than the previous value. This test is
4240	! also handled in the calling routine.
4241
4242	IMPLICIT NONE
4243
4244	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4245	ims , ime , jms , jme , kms , kme , &
4246	its , ite , jts , jte , kts , kte
4247
4248	REAL :: p_top
4249	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
4250
4251	! Local vars
4252
4253	INTEGER :: i , j , k, min_lev
4254
4255	i = its
4256	j = jts
4257	p_top = p(i,2,j)
4258	min_lev = 2
4259	DO k = 2 , kte
4260	IF ( p_top .GT. p(i,k,j) ) THEN
4261	p_top = p(i,k,j)
4262	min_lev = k
4263	END IF
4264	END DO
4265
4266	k = min_lev
4267	p_top = p(its,k,jts)
4268	DO j = jts , MIN ( jde-1 , jte )
4269	DO i = its , MIN (ide-1 , ite )
4270	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4271	p_top = MAX ( p_top , p(i,k,j) )
4272	END DO
4273	END DO
4274
4275	END SUBROUTINE find_p_top
4276
4277	!---------------------------------------------------------------------
4278
4279	SUBROUTINE t_to_theta ( t , p , p00 , &
4280	ids , ide , jds , jde , kds , kde , &
4281	ims , ime , jms , jme , kms , kme , &
4282	its , ite , jts , jte , kts , kte )
4283
4284	! Compute potential temperature from temperature and pressure.
4285
4286	IMPLICIT NONE
4287
4288	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4289	ims , ime , jms , jme , kms , kme , &
4290	its , ite , jts , jte , kts , kte
4291
4292	REAL , INTENT(IN) :: p00
4293	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
4294	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: t
4295
4296	! Local vars
4297
4298	INTEGER :: i , j , k
4299
4300	REAL , PARAMETER :: Rd = r_d
4301
4302	DO j = jts , MIN ( jde-1 , jte )
4303	DO k = kts , kte
4304	DO i = its , MIN (ide-1 , ite )
4305	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4306	t(i,k,j) = t(i,k,j) * ( p00 / p(i,k,j) ) ** (Rd / Cp)
4307	END DO
4308	END DO
4309	END DO
4310
4311	END SUBROUTINE t_to_theta
4312
4313
4314	!---------------------------------------------------------------------
4315
4316	SUBROUTINE theta_to_t ( t , p , p00 , &
4317	ids , ide , jds , jde , kds , kde , &
4318	ims , ime , jms , jme , kms , kme , &
4319	its , ite , jts , jte , kts , kte )
4320
4321	! Compute temperature from potential temp and pressure.
4322
4323	IMPLICIT NONE
4324
4325	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4326	ims , ime , jms , jme , kms , kme , &
4327	its , ite , jts , jte , kts , kte
4328
4329	REAL , INTENT(IN) :: p00
4330	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
4331	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: t
4332
4333	! Local vars
4334
4335	INTEGER :: i , j , k
4336
4337	REAL , PARAMETER :: Rd = r_d
4338	CHARACTER (LEN=80) :: mess
4339
4340	DO j = jts , MIN ( jde-1 , jte )
4341	DO k = kts , kte-1
4342	DO i = its , MIN (ide-1 , ite )
4343	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4344	if ( p(i,k,j) .NE. 0. ) then
4345	t(i,k,j) = t(i,k,j) / ( ( p00 / p(i,k,j) ) ** (Rd / Cp) )
4346	else
4347	WRITE(mess,*) 'Troubles in theta_to_t'
4348	CALL wrf_debug(0,mess)
4349	WRITE(mess,*) "i,j,k = ", i,j,k
4350	CALL wrf_debug(0,mess)
4351	WRITE(mess,*) "p(i,k,j) = ", p(i,k,j)
4352	CALL wrf_debug(0,mess)
4353	WRITE(mess,*) "t(i,k,j) = ", t(i,k,j)
4354	CALL wrf_debug(0,mess)
4355	endif
4356	END DO
4357	END DO
4358	END DO
4359
4360	END SUBROUTINE theta_to_t
4361
4362	!---------------------------------------------------------------------
4363
4364	SUBROUTINE integ_moist ( q_in , p_in , pd_out , t_in , ght_in , intq , &
4365	ids , ide , jds , jde , kds , kde , &
4366	ims , ime , jms , jme , kms , kme , &
4367	its , ite , jts , jte , kts , kte )
4368
4369	! Integrate the moisture field vertically. Mostly used to get the total
4370	! vapor pressure, which can be subtracted from the total pressure to get
4371	! the dry pressure.
4372
4373	IMPLICIT NONE
4374
4375	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4376	ims , ime , jms , jme , kms , kme , &
4377	its , ite , jts , jte , kts , kte
4378
4379	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: q_in , p_in , t_in , ght_in
4380	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pd_out
4381	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: intq
4382
4383	! Local vars
4384
4385	INTEGER :: i , j , k
4386	INTEGER , DIMENSION(ims:ime) :: level_above_sfc
4387	REAL , DIMENSION(ims:ime,jms:jme) :: psfc , tsfc , qsfc, zsfc
4388	REAL , DIMENSION(ims:ime,kms:kme) :: q , p , t , ght, pd
4389
4390	REAL :: rhobar , qbar , dz
4391	REAL :: p1 , p2 , t1 , t2 , q1 , q2 , z1, z2
4392
4393	LOGICAL :: upside_down
4394
4395	REAL , PARAMETER :: Rd = r_d
4396
4397	! Is the data upside down?
4398
4399
4400	find_valid : DO j = jts , MIN ( jde-1 , jte )
4401	DO i = its , MIN (ide-1 , ite )
4402	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4403	IF ( p_in(i,kts+1,j) .LT. p_in(i,kte,j) ) THEN
4404	upside_down = .TRUE.
4405	ELSE
4406	upside_down = .FALSE.
4407	END IF
4408	EXIT find_valid
4409	END DO
4410	END DO find_valid
4411
4412	! Get a surface value, always the first level of a 3d field.
4413
4414	DO j = jts , MIN ( jde-1 , jte )
4415	DO i = its , MIN (ide-1 , ite )
4416	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4417	psfc(i,j) = p_in(i,kts,j)
4418	tsfc(i,j) = t_in(i,kts,j)
4419	qsfc(i,j) = q_in(i,kts,j)
4420	zsfc(i,j) = ght_in(i,kts,j)
4421	END DO
4422	END DO
4423
4424	DO j = jts , MIN ( jde-1 , jte )
4425
4426	! Initialize the integrated quantity of moisture to zero.
4427
4428	DO i = its , MIN (ide-1 , ite )
4429	intq(i,j) = 0.
4430	END DO
4431
4432	IF ( upside_down ) THEN
4433	DO i = its , MIN (ide-1 , ite )
4434	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4435	p(i,kts) = p_in(i,kts,j)
4436	t(i,kts) = t_in(i,kts,j)
4437	q(i,kts) = q_in(i,kts,j)
4438	ght(i,kts) = ght_in(i,kts,j)
4439	DO k = kts+1,kte
4440	p(i,k) = p_in(i,kte+2-k,j)
4441	t(i,k) = t_in(i,kte+2-k,j)
4442	q(i,k) = q_in(i,kte+2-k,j)
4443	ght(i,k) = ght_in(i,kte+2-k,j)
4444	END DO
4445	END DO
4446	ELSE
4447	DO i = its , MIN (ide-1 , ite )
4448	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4449	DO k = kts,kte
4450	p(i,k) = p_in(i,k ,j)
4451	t(i,k) = t_in(i,k ,j)
4452	q(i,k) = q_in(i,k ,j)
4453	ght(i,k) = ght_in(i,k ,j)
4454	END DO
4455	END DO
4456	END IF
4457
4458	! Find the first level above the ground. If all of the levels are above ground, such as
4459	! a terrain following lower coordinate, then the first level above ground is index #2.
4460
4461	DO i = its , MIN (ide-1 , ite )
4462	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4463	level_above_sfc(i) = -1
4464	IF ( p(i,kts+1) .LT. psfc(i,j) ) THEN
4465	level_above_sfc(i) = kts+1
4466	ELSE
4467	find_k : DO k = kts+1,kte-1
4468	IF ( ( p(i,k )-psfc(i,j) .GE. 0. ) .AND. &
4469	( p(i,k+1)-psfc(i,j) .LT. 0. ) ) THEN
4470	level_above_sfc(i) = k+1
4471	EXIT find_k
4472	END IF
4473	END DO find_k
4474	IF ( level_above_sfc(i) .EQ. -1 ) THEN
4475	print *,'i,j = ',i,j
4476	print *,'p = ',p(i,:)
4477	print *,'p sfc = ',psfc(i,j)
4478	CALL wrf_error_fatal ( 'Could not find level above ground')
4479	END IF
4480	END IF
4481	END DO
4482
4483	DO i = its , MIN (ide-1 , ite )
4484	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4485
4486	! Account for the moisture above the ground.
4487
4488	pd(i,kte) = p(i,kte)
4489	DO k = kte-1,level_above_sfc(i),-1
4490	rhobar = ( p(i,k ) / ( Rd * t(i,k ) ) + &
4491	p(i,k+1) / ( Rd * t(i,k+1) ) ) * 0.5
4492	qbar = ( q(i,k ) + q(i,k+1) ) * 0.5
4493	dz = ght(i,k+1) - ght(i,k)
4494	intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz
4495	pd(i,k) = p(i,k) - intq(i,j)
4496	END DO
4497
4498	! Account for the moisture between the surface and the first level up.
4499
4500	IF ( ( p(i,level_above_sfc(i)-1)-psfc(i,j) .GE. 0. ) .AND. &
4501	( p(i,level_above_sfc(i) )-psfc(i,j) .LT. 0. ) .AND. &
4502	( level_above_sfc(i) .GT. kts ) ) THEN
4503	p1 = psfc(i,j)
4504	p2 = p(i,level_above_sfc(i))
4505	t1 = tsfc(i,j)
4506	t2 = t(i,level_above_sfc(i))
4507	q1 = qsfc(i,j)
4508	q2 = q(i,level_above_sfc(i))
4509	z1 = zsfc(i,j)
4510	z2 = ght(i,level_above_sfc(i))
4511	rhobar = ( p1 / ( Rd * t1 ) + &
4512	p2 / ( Rd * t2 ) ) * 0.5
4513	qbar = ( q1 + q2 ) * 0.5
4514	dz = z2 - z1
4515	IF ( dz .GT. 0.1 ) THEN
4516	intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz
4517	END IF
4518
4519	! Fix the underground values.
4520
4521	DO k = level_above_sfc(i)-1,kts+1,-1
4522	pd(i,k) = p(i,k) - intq(i,j)
4523	END DO
4524	END IF
4525	pd(i,kts) = psfc(i,j) - intq(i,j)
4526
4527	END DO
4528
4529	IF ( upside_down ) THEN
4530	DO i = its , MIN (ide-1 , ite )
4531	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4532	pd_out(i,kts,j) = pd(i,kts)
4533	DO k = kts+1,kte
4534	pd_out(i,kte+2-k,j) = pd(i,k)
4535	END DO
4536	END DO
4537	ELSE
4538	DO i = its , MIN (ide-1 , ite )
4539	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4540	DO k = kts,kte
4541	pd_out(i,k,j) = pd(i,k)
4542	END DO
4543	END DO
4544	END IF
4545
4546	END DO
4547
4548	END SUBROUTINE integ_moist
4549
4550	!---------------------------------------------------------------------
4551
4552	SUBROUTINE rh_to_mxrat2(rh, t, p, q , wrt_liquid , &
4553	qv_max_p_safe , &
4554	qv_max_flag , qv_max_value , &
4555	qv_min_p_safe , &
4556	qv_min_flag , qv_min_value , &
4557	ids , ide , jds , jde , kds , kde , &
4558	ims , ime , jms , jme , kms , kme , &
4559	its , ite , jts , jte , kts , kte )
4560
4561	! This subroutine computes mixing ratio (q, kg/kg) from basic variables
4562	! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 0-100%).
4563	! Phase transition, liquid water to ice, occurs over (0,-23) temperature range (Celcius).
4564	! Formulation used here is based on:
4565	! WMO, General meteorological standards and recommended practices,
4566	! Appendix A, WMO Technical Regulations, WMO-No. 49, corrigendum,
4567	! August 2000. --TKW 03/30/2011
4568
4569	IMPLICIT NONE
4570
4571	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4572	ims , ime , jms , jme , kms , kme , &
4573	its , ite , jts , jte , kts , kte
4574
4575	LOGICAL , INTENT(IN) :: wrt_liquid
4576
4577	REAL , INTENT(IN) :: qv_max_p_safe , qv_max_flag , qv_max_value
4578	REAL , INTENT(IN) :: qv_min_p_safe , qv_min_flag , qv_min_value
4579
4580	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
4581	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
4582	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: q
4583
4584	! Local vars
4585
4586	REAL, PARAMETER :: T0K = 273.16
4587	REAL, PARAMETER :: Tice = T0K - 23.0
4588
4589	REAL, PARAMETER :: cfe = 1.0/(23.0*23.0)
4590	REAL, PARAMETER :: eps = 0.622
4591
4592	! Coefficients for esat over liquid water
4593	REAL, PARAMETER :: cw1 = 10.79574
4594	REAL, PARAMETER :: cw2 = -5.02800
4595	REAL, PARAMETER :: cw3 = 1.50475E-4
4596	REAL, PARAMETER :: cw4 = 0.42873E-3
4597	REAL, PARAMETER :: cw5 = 0.78614
4598
4599	! Coefficients for esat over ice
4600	REAL, PARAMETER :: ci1 = -9.09685
4601	REAL, PARAMETER :: ci2 = -3.56654
4602	REAL, PARAMETER :: ci3 = 0.87682
4603	REAL, PARAMETER :: ci4 = 0.78614
4604
4605	REAL, PARAMETER :: Tn = 273.16
4606
4607	! 1 ppm is a reasonable estimate for minimum QV even for stratospheric altitudes
4608	REAL, PARAMETER :: QV_MIN = 1.e-6
4609
4610	! Maximum allowed QV is computed under the extreme condition:
4611	! Saturated at 40 degree in Celcius and 1000 hPa
4612	REAL, PARAMETER :: QV_MAX = 0.045
4613
4614	! Need to constrain WVP in the stratosphere where pressure
4615	! is low but tempearure is hot (warm)
4616	! Maximum ratio of e/p, = q/(0.622+q)
4617	REAL, PARAMETER :: EP_MAX = QV_MAX/(eps+QV_MAX)
4618
4619	INTEGER :: i , j , k
4620
4621	REAL :: ew , q1 , t1
4622	REAL :: ta, tb, pw3, pw4, pwr
4623	REAL :: es, esw, esi, wvp, pmb, wvpmax
4624
4625	DO j = jts , MIN ( jde-1 , jte )
4626	DO k = kts , kte
4627	DO i = its , MIN (ide-1 , ite )
4628	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4629	rh(i,k,j) = MIN ( MAX ( rh(i,k,j) , 0. ) , 100. )
4630	END DO
4631	END DO
4632	END DO
4633
4634	IF ( wrt_liquid ) THEN
4635	DO j = jts , MIN ( jde-1 , jte )
4636	DO k = kts , kte
4637	DO i = its , MIN (ide-1 , ite )
4638	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4639	Ta=Tn/T(i,k,j)
4640	Tb=T(i,k,j)/Tn
4641	pw3 = -8.2969*(Tb-1.0)
4642	pw4 = 4.76955*(1.0-Ta)
4643	pwr = cw1(1.0-Ta) + cw2LOG10(Tb) + cw3(1.0-10.0pw3) + cw4(10.0**pw4-1.0) + cw5
4644	es = 10.0**pwr ! Saturation WVP
4645	wvp = 0.01rh(i,k,j)es ! Actual WVP
4646	pmb = p(i,k,j)/100.
4647	wvpmax = EP_MAX*pmb ! Prevents unrealistic QV in the stratosphere
4648	wvp = MIN(wvp,wvpmax)
4649	q(i,k,j) = eps*wvp/(pmb-wvp)
4650	END DO
4651	END DO
4652	END DO
4653
4654	ELSE
4655	DO j = jts , MIN ( jde-1 , jte )
4656	DO k = kts , kte
4657	DO i = its , MIN (ide-1 , ite )
4658	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4659	Ta=Tn/T(i,k,j)
4660	Tb=T(i,k,j)/Tn
4661	IF (t(i,k,j) >= T0K) THEN ! Over liquid water
4662	pw3 = -8.2969*(Tb-1.0)
4663	pw4 = 4.76955*(1.0-Ta)
4664	pwr = cw1(1.0-Ta) + cw2LOG10(Tb) + cw3(1.0-10.0pw3) + cw4(10.0**pw4-1.0) + cw5
4665	es = 10.0**pwr
4666	wvp = 0.01rh(i,k,j)es
4667	ELSE IF (t(i,k,j) <= Tice) THEN ! Over ice
4668	pwr = ci1(Ta-1.0) + ci2LOG10(Ta) + ci3*(1.0-Tb) + ci4
4669	es = 10.0**pwr
4670	wvp = 0.01rh(i,k,j)es
4671	ELSE ! Mixed
4672	pw3 = -8.2969*(Tb-1.0)
4673	pw4 = 4.76955*(1.0-Ta)
4674	pwr = cw1(1.0-Ta) + cw2LOG10(Tb) + cw3(1.0-10.0pw3) + cw4(10.0**pw4-1.0) + cw5
4675	esw = 10.0**pwr ! Over liquid water
4676
4677	pwr = ci1(Ta-1.0) + ci2LOG10(Ta) + ci3*(1.0-Tb) + ci4
4678	esi = 10.0**pwr ! Over ice
4679
4680	es = esi + (esw-esi)cfe(T(i,k,j)-Tice)*(T(i,k,j)-Tice)
4681	wvp = 0.01rh(i,k,j)es
4682	END IF
4683	pmb = p(i,k,j)/100.
4684	wvpmax = EP_MAX*pmb ! Prevents unrealistic QV in the stratosphere
4685	wvp = MIN(wvp,wvpmax)
4686	q(i,k,j) = eps*wvp/(pmb-wvp)
4687	END DO
4688	END DO
4689	END DO
4690	END IF
4691
4692	! For pressures above a defined level, reasonable Qv values should be
4693	! a certain value or smaller. If they are larger than this, the input data
4694	! probably had "missing" RH, and we filled in some values. This is an
4695	! attempt to catch those. Also, set the minimum value for the entire
4696	! domain that is above the selected pressure level.
4697
4698	DO j = jts , MIN ( jde-1 , jte )
4699	DO k = kts , kte
4700	DO i = its , MIN (ide-1 , ite )
4701	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4702	IF ( p(i,k,j) .LT. qv_max_p_safe ) THEN
4703	IF ( q(i,k,j) .GT. qv_max_flag ) THEN
4704	q(i,k,j) = qv_max_value
4705	END IF
4706	END IF
4707	IF ( p(i,k,j) .LT. qv_min_p_safe ) THEN
4708	IF ( q(i,k,j) .LT. qv_min_flag ) THEN
4709	q(i,k,j) = qv_min_value
4710	END IF
4711	END IF
4712	END DO
4713	END DO
4714	END DO
4715
4716	END SUBROUTINE rh_to_mxrat2
4717
4718	!---------------------------------------------------------------------
4719
4720	SUBROUTINE rh_to_mxrat1(rh, t, p, q , wrt_liquid , &
4721	qv_max_p_safe , &
4722	qv_max_flag , qv_max_value , &
4723	qv_min_p_safe , &
4724	qv_min_flag , qv_min_value , &
4725	ids , ide , jds , jde , kds , kde , &
4726	ims , ime , jms , jme , kms , kme , &
4727	its , ite , jts , jte , kts , kte )
4728
4729	IMPLICIT NONE
4730
4731	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4732	ims , ime , jms , jme , kms , kme , &
4733	its , ite , jts , jte , kts , kte
4734
4735	LOGICAL , INTENT(IN) :: wrt_liquid
4736
4737	REAL , INTENT(IN) :: qv_max_p_safe , qv_max_flag , qv_max_value
4738	REAL , INTENT(IN) :: qv_min_p_safe , qv_min_flag , qv_min_value
4739
4740	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
4741	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
4742	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: q
4743
4744	! Local vars
4745
4746	INTEGER :: i , j , k
4747
4748	REAL :: ew , q1 , t1
4749
4750	REAL, PARAMETER :: T_REF = 0.0
4751	REAL, PARAMETER :: MW_AIR = 28.966
4752	REAL, PARAMETER :: MW_VAP = 18.0152
4753
4754	REAL, PARAMETER :: A0 = 6.107799961
4755	REAL, PARAMETER :: A1 = 4.436518521e-01
4756	REAL, PARAMETER :: A2 = 1.428945805e-02
4757	REAL, PARAMETER :: A3 = 2.650648471e-04
4758	REAL, PARAMETER :: A4 = 3.031240396e-06
4759	REAL, PARAMETER :: A5 = 2.034080948e-08
4760	REAL, PARAMETER :: A6 = 6.136820929e-11
4761
4762	REAL, PARAMETER :: ES0 = 6.1121
4763
4764	REAL, PARAMETER :: C1 = 9.09718
4765	REAL, PARAMETER :: C2 = 3.56654
4766	REAL, PARAMETER :: C3 = 0.876793
4767	REAL, PARAMETER :: EIS = 6.1071
4768	REAL :: RHS
4769	REAL, PARAMETER :: TF = 273.16
4770	REAL :: TK
4771
4772	REAL :: ES
4773	REAL :: QS
4774	REAL, PARAMETER :: EPS = 0.622
4775	REAL, PARAMETER :: SVP1 = 0.6112
4776	REAL, PARAMETER :: SVP2 = 17.67
4777	REAL, PARAMETER :: SVP3 = 29.65
4778	REAL, PARAMETER :: SVPT0 = 273.15
4779
4780	CHARACTER (LEN=80) :: mess
4781
4782	! This subroutine computes mixing ratio (q, kg/kg) from basic variables
4783	! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%).
4784	! The reference temperature (t_ref, C) is used to describe the temperature
4785	! at which the liquid and ice phase change occurs.
4786
4787	DO j = jts , MIN ( jde-1 , jte )
4788	DO k = kts , kte-1
4789	DO i = its , MIN (ide-1 , ite )
4790	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4791	rh(i,k,j) = MIN ( MAX ( rh(i,k,j) , 0. ) , 100. )
4792	END DO
4793	END DO
4794	END DO
4795
4796	IF ( wrt_liquid ) THEN
4797	DO j = jts , MIN ( jde-1 , jte )
4798	DO k = kts , kte-1
4799	DO i = its , MIN (ide-1 , ite )
4800	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4801
4802	! es is reduced by RH here to avoid problems in low-pressure cases
4803	if (t(i,k,j) .ne. 0.) then
4804	es=.01rh(i,k,j)svp110.EXP(svp2*(t(i,k,j)-svpt0)/(t(i,k,j)-svp3))
4805	IF (es .ge. p(i,k,j)/100.)THEN
4806	q(i,k,j)=0.0
4807	WRITE(mess,*) 'Warning: vapor pressure exceeds total pressure, setting Qv to 0'
4808	CALL wrf_debug(0,mess)
4809	ELSE
4810	q(i,k,j)=eps*es/(p(i,k,j)/100.-es)
4811	ENDIF
4812	else
4813	q(i,k,j)=0.0
4814	WRITE(mess,*) 't(i,j,k) was 0 at ', i,j,k,', setting Qv to 0'
4815	CALL wrf_debug(0,mess)
4816	endif
4817	END DO
4818	END DO
4819	END DO
4820
4821	ELSE
4822	DO j = jts , MIN ( jde-1 , jte )
4823	DO k = kts , kte-1
4824	DO i = its , MIN (ide-1 , ite )
4825	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4826
4827	t1 = t(i,k,j) - 273.16
4828
4829	! Obviously dry.
4830
4831	IF ( t1 .lt. -200. ) THEN
4832	q(i,k,j) = 0
4833
4834	ELSE
4835
4836	! First compute the ambient vapor pressure of water
4837
4838	! Liquid phase t > 0 C
4839
4840	IF ( t1 .GE. t_ref ) THEN
4841	ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6)))))
4842
4843	! Mixed phase -47 C < t < 0 C
4844
4845	ELSE IF ( ( t1 .LT. t_ref ) .AND. ( t1 .GE. -47. ) ) THEN
4846	ew = es0 * exp(17.67 * t1 / ( t1 + 243.5))
4847
4848	! Ice phase t < -47 C
4849
4850	ELSE IF ( t1 .LT. -47. ) THEN
4851	tk = t(i,k,j)
4852	rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + &
4853	c3 * (1. - tk / tf) + alog10(eis)
4854	ew = 10. ** rhs
4855
4856	END IF
4857
4858	! Now sat vap pres obtained compute local vapor pressure
4859
4860	ew = MAX ( ew , 0. ) * rh(i,k,j) * 0.01
4861
4862	! Now compute the specific humidity using the partial vapor
4863	! pressures of water vapor (ew) and dry air (p-ew). The
4864	! constants assume that the pressure is in hPa, so we divide
4865	! the pressures by 100.
4866
4867	q1 = mw_vap * ew
4868	q1 = q1 / (q1 + mw_air * (p(i,k,j)/100. - ew))
4869
4870	q(i,k,j) = q1 / (1. - q1 )
4871
4872	END IF
4873
4874	END DO
4875	END DO
4876	END DO
4877	END IF
4878
4879	! For pressures above a defined level, reasonable Qv values should be
4880	! a certain value or smaller. If they are larger than this, the input data
4881	! probably had "missing" RH, and we filled in some values. This is an
4882	! attempt to catch those. Also, set the minimum value for the entire
4883	! domain that is above the selected pressure level.
4884
4885	DO j = jts , MIN ( jde-1 , jte )
4886	DO k = kts , kte-1
4887	DO i = its , MIN (ide-1 , ite )
4888	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
4889	IF ( p(i,k,j) .LT. qv_max_p_safe ) THEN
4890	IF ( q(i,k,j) .GT. qv_max_flag ) THEN
4891	q(i,k,j) = qv_max_value
4892	END IF
4893	END IF
4894	IF ( p(i,k,j) .LT. qv_min_p_safe ) THEN
4895	IF ( q(i,k,j) .LT. qv_min_flag ) THEN
4896	q(i,k,j) = qv_min_value
4897	END IF
4898	END IF
4899	END DO
4900	END DO
4901	END DO
4902
4903	END SUBROUTINE rh_to_mxrat1
4904
4905	!---------------------------------------------------------------------
4906
4907	SUBROUTINE compute_eta ( znw , &
4908	eta_levels , max_eta , max_dz , &
4909	p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso , &
4910	ids , ide , jds , jde , kds , kde , &
4911	ims , ime , jms , jme , kms , kme , &
4912	its , ite , jts , jte , kts , kte )
4913
4914	! Compute eta levels, either using given values from the namelist (hardly
4915	! a computation, yep, I know), or assuming a constant dz above the PBL,
4916	! knowing p_top and the number of eta levels.
4917
4918	IMPLICIT NONE
4919
4920	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4921	ims , ime , jms , jme , kms , kme , &
4922	its , ite , jts , jte , kts , kte
4923	REAL , INTENT(IN) :: max_dz
4924	REAL , INTENT(IN) :: p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso
4925	INTEGER , INTENT(IN) :: max_eta
4926	REAL , DIMENSION (max_eta) , INTENT(IN) :: eta_levels
4927
4928	REAL , DIMENSION (kts:kte) , INTENT(OUT) :: znw
4929
4930	! Local vars
4931
4932	INTEGER :: k
4933	REAL :: mub , t_init , p_surf , pb, ztop, ztop_pbl , dz , temp
4934	REAL , DIMENSION(kts:kte) :: dnw
4935
4936	INTEGER , PARAMETER :: prac_levels = 17
4937	INTEGER :: loop , loop1
4938	REAL , DIMENSION(prac_levels) :: znw_prac , znu_prac , dnw_prac
4939	REAL , DIMENSION(kts:kte) :: alb , phb
4940
4941	! Gee, do the eta levels come in from the namelist?
4942
4943	IF ( ABS(eta_levels(1)+1.) .GT. 0.0000001 ) THEN
4944
4945	! Check to see if the array is oriented OK, we can easily fix an upside down oops.
4946
4947	IF ( ( ABS(eta_levels(1 )-1.) .LT. 0.0000001 ) .AND. &
4948	( ABS(eta_levels(kde)-0.) .LT. 0.0000001 ) ) THEN
4949	DO k = kds+1 , kde-1
4950	znw(k) = eta_levels(k)
4951	END DO
4952	znw( 1) = 1.
4953	znw(kde) = 0.
4954	ELSE IF ( ( ABS(eta_levels(kde)-1.) .LT. 0.0000001 ) .AND. &
4955	( ABS(eta_levels(1 )-0.) .LT. 0.0000001 ) ) THEN
4956	DO k = kds+1 , kde-1
4957	znw(k) = eta_levels(kde+1-k)
4958	END DO
4959	znw( 1) = 1.
4960	znw(kde) = 0.
4961	ELSE
4962	CALL wrf_error_fatal ( 'First eta level should be 1.0 and the last 0.0 in namelist' )
4963	END IF
4964
4965	! Check to see if the input full-level eta array is monotonic.
4966
4967	DO k = kds , kde-1
4968	IF ( znw(k) .LE. znw(k+1) ) THEN
4969	PRINT *,'eta on full levels is not monotonic'
4970	PRINT *,'eta (',k,') = ',znw(k)
4971	PRINT *,'eta (',k+1,') = ',znw(k+1)
4972	CALL wrf_error_fatal ( 'Fix non-monotonic "eta_levels" in the namelist.input file' )
4973	END IF
4974	END DO
4975
4976	! Compute eta levels assuming a constant delta z above the PBL.
4977
4978	ELSE
4979
4980	! Compute top of the atmosphere with some silly levels. We just want to
4981	! integrate to get a reasonable value for ztop. We use the planned PBL-esque
4982	! levels, and then just coarse resolution above that. We know p_top, and we
4983	! have the base state vars.
4984
4985	p_surf = p00
4986
4987	znw_prac = (/ 1.000 , 0.993 , 0.983 , 0.970 , 0.954 , 0.934 , 0.909 , &
4988	0.88 , 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /)
4989
4990	DO k = 1 , prac_levels - 1
4991	znu_prac(k) = ( znw_prac(k) + znw_prac(k+1) ) * 0.5
4992	dnw_prac(k) = znw_prac(k+1) - znw_prac(k)
4993	END DO
4994
4995	DO k = 1, prac_levels-1
4996	pb = znu_prac(k)*(p_surf - p_top) + p_top
4997	temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
4998	! temp = t00 + A*LOG(pb/p00)
4999	t_init = temp(p00/pb)*(r_d/cp) - t0
5000	alb(k) = (r_d/p1000mb)(t_init+t0)(pb/p1000mb)**cvpm
5001	END DO
5002
5003	! Base state mu is defined as base state surface pressure minus p_top
5004
5005	mub = p_surf - p_top
5006
5007	! Integrate base geopotential, starting at terrain elevation.
5008
5009	phb(1) = 0.
5010	DO k = 2,prac_levels
5011	phb(k) = phb(k-1) - dnw_prac(k-1)mubalb(k-1)
5012	END DO
5013
5014	! So, now we know the model top in meters. Get the average depth above the PBL
5015	! of each of the remaining levels. We are going for a constant delta z thickness.
5016
5017	ztop = phb(prac_levels) / g
5018	ztop_pbl = phb(8 ) / g
5019	dz = ( ztop - ztop_pbl ) / REAL ( kde - 8 )
5020
5021	! Standard levels near the surface so no one gets in trouble.
5022
5023	DO k = 1 , 8
5024	znw(k) = znw_prac(k)
5025	END DO
5026
5027	! Using d phb(k)/ d eta(k) = -mub * alb(k), eqn 2.9
5028	! Skamarock et al, NCAR TN 468. Use full levels, so
5029	! use twice the thickness.
5030
5031	DO k = 8, kte-1-2
5032	pb = znw(k) * (p_surf - p_top) + p_top
5033	temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
5034	! temp = t00 + A*LOG(pb/p00)
5035	t_init = temp(p00/pb)*(r_d/cp) - t0
5036	alb(k) = (r_d/p1000mb)(t_init+t0)(pb/p1000mb)**cvpm
5037	znw(k+1) = znw(k) - dzg / ( mubalb(k) )
5038	END DO
5039	znw(kte-2) = 0.000
5040
5041	! There is some iteration. We want the top level, ztop, to be
5042	! consistent with the delta z, and we want the half level values
5043	! to be consistent with the eta levels. The inner loop to 10 gets
5044	! the eta levels very accurately, but has a residual at the top, due
5045	! to dz changing. We reset dz five times, and then things seem OK.
5046
5047	DO loop1 = 1 , 5
5048	DO loop = 1 , 10
5049	DO k = 8, kte-1-2
5050	pb = (znw(k)+znw(k+1))0.5 (p_surf - p_top) + p_top
5051	temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
5052	! temp = t00 + A*LOG(pb/p00)
5053	t_init = temp(p00/pb)*(r_d/cp) - t0
5054	alb(k) = (r_d/p1000mb)(t_init+t0)(pb/p1000mb)**cvpm
5055	znw(k+1) = znw(k) - dzg / ( mubalb(k) )
5056	END DO
5057	IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN
5058	print *,'Converged znw(kte) should be about 0.0 = ',znw(kte-2)
5059	END IF
5060	znw(kte-2) = 0.000
5061	END DO
5062
5063	! Here is where we check the eta levels values we just computed.
5064
5065	DO k = 1, kde-1-2
5066	pb = (znw(k)+znw(k+1))0.5 (p_surf - p_top) + p_top
5067	temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
5068	! temp = t00 + A*LOG(pb/p00)
5069	t_init = temp(p00/pb)*(r_d/cp) - t0
5070	alb(k) = (r_d/p1000mb)(t_init+t0)(pb/p1000mb)**cvpm
5071	END DO
5072
5073	phb(1) = 0.
5074	DO k = 2,kde-2
5075	phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1)
5076	END DO
5077
5078	! Reset the model top and the dz, and iterate.
5079
5080	ztop = phb(kde-2)/g
5081	ztop_pbl = phb(8)/g
5082	dz = ( ztop - ztop_pbl ) / REAL ( (kde-2) - 8 )
5083	END DO
5084
5085	IF ( dz .GT. max_dz ) THEN
5086	print *,'z (m) = ',phb(1)/g
5087	do k = 2 ,kte-2
5088	print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g
5089	end do
5090	print *,'dz (m) above fixed eta levels = ',dz
5091	print *,'namelist max_dz (m) = ',max_dz
5092	print *,'namelist p_top (Pa) = ',p_top
5093	CALL wrf_debug ( 0, 'You need one of three things:' )
5094	CALL wrf_debug ( 0, '1) More eta levels to reduce the dz: e_vert' )
5095	CALL wrf_debug ( 0, '2) A lower p_top so your total height is reduced: p_top_requested')
5096	CALL wrf_debug ( 0, '3) Increase the maximum allowable eta thickness: max_dz')
5097	CALL wrf_debug ( 0, 'All are namelist options')
5098	CALL wrf_error_fatal ( 'dz above fixed eta levels is too large')
5099	END IF
5100
5101	! Add those 2 levels back into the middle, just above the 8 levels
5102	! that semi define a boundary layer. After we open up the levels,
5103	! then we just linearly interpolate in znw. So now levels 1-8 are
5104	! specified as the fixed boundary layer levels given in this routine.
5105	! The top levels, 12 through kte are those computed. The middle
5106	! levels 9, 10, and 11 are equi-spaced in znw, and are each 1/2 the
5107	! the znw thickness of levels 11 through 12.
5108
5109	DO k = kte-2 , 9 , -1
5110	znw(k+2) = znw(k)
5111	END DO
5112
5113	znw( 9) = 0.75 * znw( 8) + 0.25 * znw(12)
5114	znw(10) = 0.50 * znw( 8) + 0.50 * znw(12)
5115	znw(11) = 0.25 * znw( 8) + 0.75 * znw(12)
5116
5117	END IF
5118
5119	END SUBROUTINE compute_eta
5120
5121	!---------------------------------------------------------------------
5122
5123	SUBROUTINE monthly_min_max ( field_in , field_min , field_max , &
5124	ids , ide , jds , jde , kds , kde , &
5125	ims , ime , jms , jme , kms , kme , &
5126	its , ite , jts , jte , kts , kte )
5127
5128	! Plow through each month, find the max, min values for each i,j.
5129
5130	IMPLICIT NONE
5131
5132	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5133	ims , ime , jms , jme , kms , kme , &
5134	its , ite , jts , jte , kts , kte
5135
5136	REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in
5137	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max
5138
5139	! Local vars
5140
5141	INTEGER :: i , j , l
5142	REAL :: minner , maxxer
5143
5144	DO j = jts , MIN(jde-1,jte)
5145	DO i = its , MIN(ide-1,ite)
5146	minner = field_in(i,1,j)
5147	maxxer = field_in(i,1,j)
5148	DO l = 2 , 12
5149	IF ( field_in(i,l,j) .LT. minner ) THEN
5150	minner = field_in(i,l,j)
5151	END IF
5152	IF ( field_in(i,l,j) .GT. maxxer ) THEN
5153	maxxer = field_in(i,l,j)
5154	END IF
5155	END DO
5156	field_min(i,j) = minner
5157	field_max(i,j) = maxxer
5158	END DO
5159	END DO
5160
5161	END SUBROUTINE monthly_min_max
5162
5163	!---------------------------------------------------------------------
5164
5165	SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , &
5166	ids , ide , jds , jde , kds , kde , &
5167	ims , ime , jms , jme , kms , kme , &
5168	its , ite , jts , jte , kts , kte )
5169
5170	! Linrarly in time interpolate data to a current valid time. The data is
5171	! assumed to come in "monthly", valid at the 15th of every month.
5172
5173	IMPLICIT NONE
5174
5175	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5176	ims , ime , jms , jme , kms , kme , &
5177	its , ite , jts , jte , kts , kte
5178
5179	CHARACTER (LEN=24) , INTENT(IN) :: date_str
5180	REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in
5181	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out
5182
5183	! Local vars
5184
5185	INTEGER :: i , j , l
5186	INTEGER , DIMENSION(0:13) :: middle
5187	INTEGER :: target_julyr , target_julday , target_date
5188	INTEGER :: julyr , julday , int_month , month1 , month2
5189	REAL :: gmt
5190	CHARACTER (LEN=4) :: yr
5191	CHARACTER (LEN=2) :: mon , day15
5192
5193
5194	WRITE(day15,FMT='(I2.2)') 15
5195	DO l = 1 , 12
5196	WRITE(mon,FMT='(I2.2)') l
5197	CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt )
5198	middle(l) = julyr*1000 + julday
5199	END DO
5200
5201	l = 0
5202	middle(l) = middle( 1) - 31
5203
5204	l = 13
5205	middle(l) = middle(12) + 31
5206
5207	CALL get_julgmt ( date_str , target_julyr , target_julday , gmt )
5208	target_date = target_julyr * 1000 + target_julday
5209	find_month : DO l = 0 , 12
5210	IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN
5211	DO j = jts , MIN ( jde-1 , jte )
5212	DO i = its , MIN (ide-1 , ite )
5213	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5214	int_month = l
5215	IF ( ( int_month .EQ. 0 ) .OR. ( int_month .EQ. 12 ) ) THEN
5216	month1 = 12
5217	month2 = 1
5218	ELSE
5219	month1 = int_month
5220	month2 = month1 + 1
5221	END IF
5222	field_out(i,j) = ( field_in(i,month2,j) * ( target_date - middle(l) ) + &
5223	field_in(i,month1,j) * ( middle(l+1) - target_date ) ) / &
5224	( middle(l+1) - middle(l) )
5225	END DO
5226	END DO
5227	EXIT find_month
5228	END IF
5229	END DO find_month
5230
5231	END SUBROUTINE monthly_interp_to_date
5232
5233	!---------------------------------------------------------------------
5234
5235	SUBROUTINE sfcprs (t, q, height, pslv, ter, avgsfct, p, &
5236	psfc, ez_method, &
5237	ids , ide , jds , jde , kds , kde , &
5238	ims , ime , jms , jme , kms , kme , &
5239	its , ite , jts , jte , kts , kte )
5240
5241
5242	! Computes the surface pressure using the input height,
5243	! temperature and q (already computed from relative
5244	! humidity) on p surfaces. Sea level pressure is used
5245	! to extrapolate a first guess.
5246
5247	IMPLICIT NONE
5248
5249	REAL, PARAMETER :: gamma = 6.5E-3
5250	REAL, PARAMETER :: pconst = 10000.0
5251	REAL, PARAMETER :: Rd = r_d
5252	REAL, PARAMETER :: TC = svpt0 + 17.5
5253
5254	REAL, PARAMETER :: gammarg = gamma * Rd / g
5255	REAL, PARAMETER :: rov2 = Rd / 2.
5256
5257	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5258	ims , ime , jms , jme , kms , kme , &
5259	its , ite , jts , jte , kts , kte
5260	LOGICAL , INTENT ( IN ) :: ez_method
5261
5262	REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p
5263	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: pslv , ter, avgsfct
5264	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
5265
5266	INTEGER :: i
5267	INTEGER :: j
5268	INTEGER :: k
5269	INTEGER , DIMENSION (its:ite,jts:jte) :: k500 , k700 , k850
5270
5271	LOGICAL :: l1
5272	LOGICAL :: l2
5273	LOGICAL :: l3
5274	LOGICAL :: OK
5275
5276	REAL :: gamma78 ( its:ite,jts:jte )
5277	REAL :: gamma57 ( its:ite,jts:jte )
5278	REAL :: ht ( its:ite,jts:jte )
5279	REAL :: p1 ( its:ite,jts:jte )
5280	REAL :: t1 ( its:ite,jts:jte )
5281	REAL :: t500 ( its:ite,jts:jte )
5282	REAL :: t700 ( its:ite,jts:jte )
5283	REAL :: t850 ( its:ite,jts:jte )
5284	REAL :: tfixed ( its:ite,jts:jte )
5285	REAL :: tsfc ( its:ite,jts:jte )
5286	REAL :: tslv ( its:ite,jts:jte )
5287
5288	! We either compute the surface pressure from a time averaged surface temperature
5289	! (what we will call the "easy way"), or we try to remove the diurnal impact on the
5290	! surface temperature (what we will call the "other way"). Both are essentially
5291	! corrections to a sea level pressure with a high-resolution topography field.
5292
5293	IF ( ez_method ) THEN
5294
5295	DO j = jts , MIN(jde-1,jte)
5296	DO i = its , MIN(ide-1,ite)
5297	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5298	psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / avgsfct(i,j) ) ** ( - g / ( Rd * gamma ) )
5299	END DO
5300	END DO
5301
5302	ELSE
5303
5304	! Find the locations of the 850, 700 and 500 mb levels.
5305
5306	k850 = 0 ! find k at: P=850
5307	k700 = 0 ! P=700
5308	k500 = 0 ! P=500
5309
5310	i = its
5311	j = jts
5312	DO k = kts+1 , kte
5313	IF (NINT(p(i,k,j)) .EQ. 85000) THEN
5314	k850(i,j) = k
5315	ELSE IF (NINT(p(i,k,j)) .EQ. 70000) THEN
5316	k700(i,j) = k
5317	ELSE IF (NINT(p(i,k,j)) .EQ. 50000) THEN
5318	k500(i,j) = k
5319	END IF
5320	END DO
5321
5322	IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN
5323
5324	DO j = jts , MIN(jde-1,jte)
5325	DO i = its , MIN(ide-1,ite)
5326	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5327	psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / t(i,1,j) ) ** ( - g / ( Rd * gamma ) )
5328	END DO
5329	END DO
5330
5331	RETURN
5332	#if 0
5333
5334	! Possibly it is just that we have a generalized vertical coord, so we do not
5335	! have the values exactly. Do a simple assignment to a close vertical level.
5336
5337	DO j = jts , MIN(jde-1,jte)
5338	DO i = its , MIN(ide-1,ite)
5339	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5340	DO k = kts+1 , kte-1
5341	IF ( ( p(i,k,j) - 85000. ) * ( p(i,k+1,j) - 85000. ) .LE. 0.0 ) THEN
5342	k850(i,j) = k
5343	END IF
5344	IF ( ( p(i,k,j) - 70000. ) * ( p(i,k+1,j) - 70000. ) .LE. 0.0 ) THEN
5345	k700(i,j) = k
5346	END IF
5347	IF ( ( p(i,k,j) - 50000. ) * ( p(i,k+1,j) - 50000. ) .LE. 0.0 ) THEN
5348	k500(i,j) = k
5349	END IF
5350	END DO
5351	END DO
5352	END DO
5353
5354	! If we still do not have the k levels, punt. I mean, we did try.
5355
5356	OK = .TRUE.
5357	DO j = jts , MIN(jde-1,jte)
5358	DO i = its , MIN(ide-1,ite)
5359	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5360	IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN
5361	OK = .FALSE.
5362	PRINT '(A)','(i,j) = ',i,j,' Error in finding p level for 850, 700 or 500 hPa.'
5363	DO K = kts+1 , kte
5364	PRINT '(A,I3,A,F10.2,A)','K = ',k,' PRESSURE = ',p(i,k,j),' Pa'
5365	END DO
5366	PRINT '(A)','Expected 850, 700, and 500 mb values, at least.'
5367	END IF
5368	END DO
5369	END DO
5370	IF ( .NOT. OK ) THEN
5371	CALL wrf_error_fatal ( 'wrong pressure levels' )
5372	END IF
5373	#endif
5374
5375	! We are here if the data is isobaric and we found the levels for 850, 700,
5376	! and 500 mb right off the bat.
5377
5378	ELSE
5379	DO j = jts , MIN(jde-1,jte)
5380	DO i = its , MIN(ide-1,ite)
5381	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5382	k850(i,j) = k850(its,jts)
5383	k700(i,j) = k700(its,jts)
5384	k500(i,j) = k500(its,jts)
5385	END DO
5386	END DO
5387	END IF
5388
5389	! The 850 hPa level of geopotential height is called something special.
5390
5391	DO j = jts , MIN(jde-1,jte)
5392	DO i = its , MIN(ide-1,ite)
5393	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5394	ht(i,j) = height(i,k850(i,j),j)
5395	END DO
5396	END DO
5397
5398	! The variable ht is now -ter/ht(850 hPa). The plot thickens.
5399
5400	DO j = jts , MIN(jde-1,jte)
5401	DO i = its , MIN(ide-1,ite)
5402	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5403	ht(i,j) = -ter(i,j) / ht(i,j)
5404	END DO
5405	END DO
5406
5407	! Make an isothermal assumption to get a first guess at the surface
5408	! pressure. This is to tell us which levels to use for the lapse
5409	! rates in a bit.
5410
5411	DO j = jts , MIN(jde-1,jte)
5412	DO i = its , MIN(ide-1,ite)
5413	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5414	psfc(i,j) = pslv(i,j) * (pslv(i,j) / p(i,k850(i,j),j)) ** ht(i,j)
5415	END DO
5416	END DO
5417
5418	! Get a pressure more than pconst Pa above the surface - p1. The
5419	! p1 is the top of the level that we will use for our lapse rate
5420	! computations.
5421
5422	DO j = jts , MIN(jde-1,jte)
5423	DO i = its , MIN(ide-1,ite)
5424	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5425	IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN
5426	p1(i,j) = 85000.
5427	ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0. ) THEN
5428	p1(i,j) = psfc(i,j) - pconst
5429	ELSE
5430	p1(i,j) = 50000.
5431	END IF
5432	END DO
5433	END DO
5434
5435	! Compute virtual temperatures for k850, k700, and k500 layers. Now
5436	! you see why we wanted Q on pressure levels, it all is beginning
5437	! to make sense.
5438
5439	DO j = jts , MIN(jde-1,jte)
5440	DO i = its , MIN(ide-1,ite)
5441	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5442	t850(i,j) = t(i,k850(i,j),j) * (1. + 0.608 * q(i,k850(i,j),j))
5443	t700(i,j) = t(i,k700(i,j),j) * (1. + 0.608 * q(i,k700(i,j),j))
5444	t500(i,j) = t(i,k500(i,j),j) * (1. + 0.608 * q(i,k500(i,j),j))
5445	END DO
5446	END DO
5447
5448	! Compute lapse rates between these three levels. These are
5449	! environmental values for each (i,j).
5450
5451	DO j = jts , MIN(jde-1,jte)
5452	DO i = its , MIN(ide-1,ite)
5453	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5454	gamma78(i,j) = ALOG(t850(i,j) / t700(i,j)) / ALOG (p(i,k850(i,j),j) / p(i,k700(i,j),j) )
5455	gamma57(i,j) = ALOG(t700(i,j) / t500(i,j)) / ALOG (p(i,k700(i,j),j) / p(i,k500(i,j),j) )
5456	END DO
5457	END DO
5458
5459	DO j = jts , MIN(jde-1,jte)
5460	DO i = its , MIN(ide-1,ite)
5461	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5462	IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN
5463	t1(i,j) = t850(i,j)
5464	ELSE IF ( ( psfc(i,j) - 85000. ) .GE. 0. ) THEN
5465	t1(i,j) = t700(i,j) * (p1(i,j) / (p(i,k700(i,j),j))) ** gamma78(i,j)
5466	ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0.) THEN
5467	t1(i,j) = t500(i,j) * (p1(i,j) / (p(i,k500(i,j),j))) ** gamma57(i,j)
5468	ELSE
5469	t1(i,j) = t500(i,j)
5470	ENDIF
5471	END DO
5472	END DO
5473
5474	! From our temperature way up in the air, we extrapolate down to
5475	! the sea level to get a guess at the sea level temperature.
5476
5477	DO j = jts , MIN(jde-1,jte)
5478	DO i = its , MIN(ide-1,ite)
5479	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5480	tslv(i,j) = t1(i,j) * (pslv(i,j) / p1(i,j)) ** gammarg
5481	END DO
5482	END DO
5483
5484	! The new surface temperature is computed from the with new sea level
5485	! temperature, just using the elevation and a lapse rate. This lapse
5486	! rate is -6.5 K/km.
5487
5488	DO j = jts , MIN(jde-1,jte)
5489	DO i = its , MIN(ide-1,ite)
5490	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5491	tsfc(i,j) = tslv(i,j) - gamma * ter(i,j)
5492	END DO
5493	END DO
5494
5495	! A small correction to the sea-level temperature, in case it is too warm.
5496
5497	DO j = jts , MIN(jde-1,jte)
5498	DO i = its , MIN(ide-1,ite)
5499	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5500	tfixed(i,j) = tc - 0.005 * (tsfc(i,j) - tc) ** 2
5501	END DO
5502	END DO
5503
5504	DO j = jts , MIN(jde-1,jte)
5505	DO i = its , MIN(ide-1,ite)
5506	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5507	l1 = tslv(i,j) .LT. tc
5508	l2 = tsfc(i,j) .LE. tc
5509	l3 = .NOT. l1
5510	IF ( l2 .AND. l3 ) THEN
5511	tslv(i,j) = tc
5512	ELSE IF ( ( .NOT. l2 ) .AND. l3 ) THEN
5513	tslv(i,j) = tfixed(i,j)
5514	END IF
5515	END DO
5516	END DO
5517
5518	! Finally, we can get to the surface pressure.
5519
5520	DO j = jts , MIN(jde-1,jte)
5521	DO i = its , MIN(ide-1,ite)
5522	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5523	p1(i,j) = - ter(i,j) * g / ( rov2 * ( tsfc(i,j) + tslv(i,j) ) )
5524	psfc(i,j) = pslv(i,j) * EXP ( p1(i,j) )
5525	END DO
5526	END DO
5527
5528	END IF
5529
5530	! Surface pressure and sea-level pressure are the same at sea level.
5531
5532	! DO j = jts , MIN(jde-1,jte)
5533	! DO i = its , MIN(ide-1,ite)
5534	! IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5535	! IF ( ABS ( ter(i,j) ) .LT. 0.1 ) THEN
5536	! psfc(i,j) = pslv(i,j)
5537	! END IF
5538	! END DO
5539	! END DO
5540
5541	END SUBROUTINE sfcprs
5542
5543	!---------------------------------------------------------------------
5544
5545	SUBROUTINE sfcprs2(t, q, height, psfc_in, ter, avgsfct, p, &
5546	psfc, ez_method, &
5547	ids , ide , jds , jde , kds , kde , &
5548	ims , ime , jms , jme , kms , kme , &
5549	its , ite , jts , jte , kts , kte )
5550
5551
5552	! Computes the surface pressure using the input height,
5553	! temperature and q (already computed from relative
5554	! humidity) on p surfaces. Sea level pressure is used
5555	! to extrapolate a first guess.
5556
5557	IMPLICIT NONE
5558
5559	REAL, PARAMETER :: Rd = r_d
5560
5561	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5562	ims , ime , jms , jme , kms , kme , &
5563	its , ite , jts , jte , kts , kte
5564	LOGICAL , INTENT ( IN ) :: ez_method
5565
5566	REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p
5567	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: psfc_in , ter, avgsfct
5568	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
5569
5570	INTEGER :: i
5571	INTEGER :: j
5572	INTEGER :: k
5573
5574	REAL :: tv_sfc_avg , tv_sfc , del_z
5575
5576	! Compute the new surface pressure from the old surface pressure, and a
5577	! known change in elevation at the surface.
5578
5579	! del_z = diff in surface topo, lo-res vs hi-res
5580	! psfc = psfc_in * exp ( g del_z / (Rd Tv_sfc ) )
5581
5582
5583	IF ( ez_method ) THEN
5584	DO j = jts , MIN(jde-1,jte)
5585	DO i = its , MIN(ide-1,ite)
5586	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5587	tv_sfc_avg = avgsfct(i,j) * (1. + 0.608 * q(i,1,j))
5588	del_z = height(i,1,j) - ter(i,j)
5589	psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc_avg ) )
5590	END DO
5591	END DO
5592	ELSE
5593	DO j = jts , MIN(jde-1,jte)
5594	DO i = its , MIN(ide-1,ite)
5595	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5596	tv_sfc = t(i,1,j) * (1. + 0.608 * q(i,1,j))
5597	del_z = height(i,1,j) - ter(i,j)
5598	psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc ) )
5599	END DO
5600	END DO
5601	END IF
5602
5603	END SUBROUTINE sfcprs2
5604
5605	!---------------------------------------------------------------------
5606
5607	SUBROUTINE sfcprs3( height , p , ter , slp , psfc , &
5608	ids , ide , jds , jde , kds , kde , &
5609	ims , ime , jms , jme , kms , kme , &
5610	its , ite , jts , jte , kts , kte )
5611
5612	! Computes the surface pressure by vertically interpolating
5613	! linearly (or log) in z the pressure, to the targeted topography.
5614
5615	IMPLICIT NONE
5616
5617	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5618	ims , ime , jms , jme , kms , kme , &
5619	its , ite , jts , jte , kts , kte
5620
5621	REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: height, p
5622	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: ter , slp
5623	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
5624
5625	INTEGER :: i
5626	INTEGER :: j
5627	INTEGER :: k
5628
5629	LOGICAL :: found_loc
5630
5631	REAL :: zl , zu , pl , pu , zm
5632
5633	! Loop over each grid point
5634
5635	DO j = jts , MIN(jde-1,jte)
5636	DO i = its , MIN(ide-1,ite)
5637	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
5638
5639	! Special case where near the ocean level. Assume that the SLP is a good value.
5640
5641	IF ( ter(i,j) .LT. 50 ) THEN
5642	psfc(i,j) = slp(i,j) + ( p(i,2,j)-p(i,3,j) ) / ( height(i,2,j)-height(i,3,j) ) * ter(i,j)
5643	CYCLE
5644	END IF
5645
5646	! Find the trapping levels
5647
5648	found_loc = .FALSE.
5649
5650	! Normal sort of scenario - the model topography is somewhere between
5651	! the height values of 1000 mb and the top of the model.
5652
5653	found_k_loc : DO k = kts+1 , kte-2
5654	IF ( ( height(i,k ,j) .LE. ter(i,j) ) .AND. &
5655	( height(i,k+1,j) .GT. ter(i,j) ) ) THEN
5656	zl = height(i,k ,j)
5657	zu = height(i,k+1,j)
5658	zm = ter(i,j)
5659	pl = p(i,k ,j)
5660	pu = p(i,k+1,j)
5661	psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
5662	found_loc = .TRUE.
5663	EXIT found_k_loc
5664	END IF
5665	END DO found_k_loc
5666
5667	! Interpolate betwixt slp and the first isobaric level above - this is probably the
5668	! usual thing over the ocean.
5669
5670	IF ( .NOT. found_loc ) THEN
5671	IF ( slp(i,j) .GE. p(i,2,j) ) THEN
5672	zl = 0.
5673	zu = height(i,3,j)
5674	zm = ter(i,j)
5675	pl = slp(i,j)
5676	pu = p(i,3,j)
5677	psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
5678	found_loc = .TRUE.
5679	ELSE
5680	found_slp_loc : DO k = kts+1 , kte-3
5681	IF ( ( slp(i,j) .GE. p(i,k+1,j) ) .AND. &
5682	( slp(i,j) .LT. p(i,k ,j) ) ) THEN
5683	zl = 0.
5684	zu = height(i,k+1,j)
5685	zm = ter(i,j)
5686	pl = slp(i,j)
5687	pu = p(i,k+1,j)
5688	psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
5689	found_loc = .TRUE.
5690	EXIT found_slp_loc
5691	END IF
5692	END DO found_slp_loc
5693	END IF
5694	END IF
5695
5696	! Did we do what we wanted done.
5697
5698	IF ( .NOT. found_loc ) THEN
5699	print *,'i,j = ',i,j
5700	print *,'p column = ',p(i,2:,j)
5701	print *,'z column = ',height(i,2:,j)
5702	print *,'model topo = ',ter(i,j)
5703	CALL wrf_error_fatal ( ' probs with sfc p computation ' )
5704	END IF
5705
5706	END DO
5707	END DO
5708
5709	END SUBROUTINE sfcprs3
5710	!---------------------------------------------------------------------
5711
5712	SUBROUTINE filter_topo ( ht_in , xlat , msftx , fft_filter_lat , &
5713	ids , ide , jds , jde , kds , kde , &
5714	ims , ime , jms , jme , kms , kme , &
5715	its , ite , jts , jte , kts , kte )
5716
5717	IMPLICIT NONE
5718
5719	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5720	ims , ime , jms , jme , kms , kme , &
5721	its , ite , jts , jte , kts , kte
5722
5723	REAL , INTENT(IN) :: fft_filter_lat
5724	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: ht_in
5725	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: xlat , msftx
5726
5727
5728	! Local vars
5729
5730	INTEGER :: i , j , j_lat_pos , j_lat_neg
5731	INTEGER :: i_kicker , ik , i1, i2, i3, i4
5732	REAL :: length_scale , sum
5733	REAL , DIMENSION(its:ite,jts:jte) :: ht_out
5734
5735	! The filtering is a simple average on a latitude loop. Possibly a LONG list of
5736	! numbers. We assume that ALL of the 2d arrays have been transposed so that
5737	! each patch has the entire domain size of the i-dim local.
5738
5739	IF ( ( its .NE. ids ) .OR. ( ite .NE. ide ) ) THEN
5740	CALL wrf_error_fatal ( 'filtering assumes all values on X' )
5741	END IF
5742
5743	! Starting at the south pole, we find where the
5744	! grid distance is big enough, then go back a point. Continuing to the
5745	! north pole, we find the first small grid distance. These are the
5746	! computational latitude loops and the associated computational poles.
5747
5748	j_lat_neg = 0
5749	j_lat_pos = jde + 1
5750	loop_neg : DO j = jts , MIN(jde-1,jte)
5751	IF ( xlat(its,j) .LT. 0.0 ) THEN
5752	IF ( ABS(xlat(its,j)) .LT. fft_filter_lat ) THEN
5753	j_lat_neg = j - 1
5754	EXIT loop_neg
5755	END IF
5756	END IF
5757	END DO loop_neg
5758
5759	loop_pos : DO j = jts , MIN(jde-1,jte)
5760	IF ( xlat(its,j) .GT. 0.0 ) THEN
5761	IF ( xlat(its,j) .GE. fft_filter_lat ) THEN
5762	j_lat_pos = j
5763	EXIT loop_pos
5764	END IF
5765	END IF
5766	END DO loop_pos
5767
5768	! Set output values to initial input topo values for whole patch.
5769
5770	DO j = jts , MIN(jde-1,jte)
5771	DO i = its , MIN(ide-1,ite)
5772	ht_out(i,j) = ht_in(i,j)
5773	END DO
5774	END DO
5775
5776	! Filter the topo at the negative lats.
5777
5778	DO j = j_lat_neg , jts , -1
5779	i_kicker = MIN( MAX ( NINT(msftx(its,j)) , 1 ) , (ide - ids) / 2 )
5780	print *,'j = ' , j, ', kicker = ',i_kicker
5781	DO i = its , MIN(ide-1,ite)
5782	IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
5783	sum = 0.0
5784	DO ik = 1 , i_kicker
5785	sum = sum + ht_in(i+ik,j) + ht_in(i-ik,j)
5786	END DO
5787	ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
5788	ELSE IF ( ( i - i_kicker .LT. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
5789	sum = 0.0
5790	DO ik = 1 , i_kicker
5791	sum = sum + ht_in(i+ik,j)
5792	END DO
5793	i1 = i - i_kicker + ide -1
5794	i2 = ide-1
5795	i3 = ids
5796	i4 = i-1
5797	DO ik = i1 , i2
5798	sum = sum + ht_in(ik,j)
5799	END DO
5800	DO ik = i3 , i4
5801	sum = sum + ht_in(ik,j)
5802	END DO
5803	ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
5804	ELSE IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .GT. ide-1 ) ) THEN
5805	sum = 0.0
5806	DO ik = 1 , i_kicker
5807	sum = sum + ht_in(i-ik,j)
5808	END DO
5809	i1 = i+1
5810	i2 = ide-1
5811	i3 = ids
5812	i4 = ids + ( i_kicker+i ) - ide
5813	DO ik = i1 , i2
5814	sum = sum + ht_in(ik,j)
5815	END DO
5816	DO ik = i3 , i4
5817	sum = sum + ht_in(ik,j)
5818	END DO
5819	ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
5820	END IF
5821	END DO
5822	END DO
5823
5824	! Filter the topo at the positive lats.
5825
5826	DO j = j_lat_pos , MIN(jde-1,jte)
5827	i_kicker = MIN( MAX ( NINT(msftx(its,j)) , 1 ) , (ide - ids) / 2 )
5828	print *,'j = ' , j, ', kicker = ',i_kicker
5829	DO i = its , MIN(ide-1,ite)
5830	IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
5831	sum = 0.0
5832	DO ik = 1 , i_kicker
5833	sum = sum + ht_in(i+ik,j) + ht_in(i-ik,j)
5834	END DO
5835	ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
5836	ELSE IF ( ( i - i_kicker .LT. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
5837	sum = 0.0
5838	DO ik = 1 , i_kicker
5839	sum = sum + ht_in(i+ik,j)
5840	END DO
5841	i1 = i - i_kicker + ide -1
5842	i2 = ide-1
5843	i3 = ids
5844	i4 = i-1
5845	DO ik = i1 , i2
5846	sum = sum + ht_in(ik,j)
5847	END DO
5848	DO ik = i3 , i4
5849	sum = sum + ht_in(ik,j)
5850	END DO
5851	ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
5852	ELSE IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .GT. ide-1 ) ) THEN
5853	sum = 0.0
5854	DO ik = 1 , i_kicker
5855	sum = sum + ht_in(i-ik,j)
5856	END DO
5857	i1 = i+1
5858	i2 = ide-1
5859	i3 = ids
5860	i4 = ids + ( i_kicker+i ) - ide
5861	DO ik = i1 , i2
5862	sum = sum + ht_in(ik,j)
5863	END DO
5864	DO ik = i3 , i4
5865	sum = sum + ht_in(ik,j)
5866	END DO
5867	ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
5868	END IF
5869	END DO
5870	END DO
5871
5872	! Set output values to initial input topo values for whole patch.
5873
5874	DO j = jts , MIN(jde-1,jte)
5875	DO i = its , MIN(ide-1,ite)
5876	ht_in(i,j) = ht_out(i,j)
5877	END DO
5878	END DO
5879
5880	END SUBROUTINE filter_topo
5881
5882	!---------------------------------------------------------------------
5883
5884	SUBROUTINE init_module_initialize
5885	END SUBROUTINE init_module_initialize
5886
5887	!---------------------------------------------------------------------
5888
5889	END MODULE module_initialize_real
5890	#endif

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

Download in other formats:

Original Format