Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_initialize_real.F @ 3709

Last change on this file since 3709 was 2759, checked in by aslmd, 3 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 193.2 KB

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

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

Download in other formats:

Original Format