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