module_initialize_real.F

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