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