| blob | 4ea1bb6377255b85d63d2960d9f79b825762869e |
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 , tiso
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 , tiso )
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 IF ( flag_tavgsfc .EQ. 1 ) THEN
252 we_have_tavgsfc = .TRUE.
253 ELSE
254 we_have_tavgsfc = .FALSE.
255 END IF
257 ! Is there any vertical interpolation to do? The "old" data comes in on the correct
258 ! vertical locations already.
260 IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
262 ! Variables that are named differently between SI and WPS.
264 DO j = jts, MIN(jte,jde-1)
265 DO i = its, MIN(ite,ide-1)
266 grid%tsk(i,j) = grid%tsk_gc(i,j)
267 grid%tmn(i,j) = grid%tmn_gc(i,j)
268 grid%xlat(i,j) = grid%xlat_gc(i,j)
269 grid%xlong(i,j) = grid%xlong_gc(i,j)
270 grid%ht(i,j) = grid%ht_gc(i,j)
271 END DO
272 END DO
274 ! A user could request that the most coarse grid has the
275 ! topography along the outer boundary smoothed. This smoothing
276 ! is similar to the coarse/nest interface. The outer rows and
277 ! cols come from the existing large scale topo, and then the
278 ! next several rows/cols are a linear ramp of the large scale
279 ! model and the hi-res topo from WPS. We only do this for the
280 ! coarse grid since we are going to make the interface consistent
281 ! in the model betwixt the CG and FG domains.
283 IF ( ( config_flags%smooth_cg_topo ) .AND. &
284 ( grid%id .EQ. 1 ) .AND. &
285 ( flag_soilhgt .EQ. 1) ) THEN
286 CALL blend_terrain ( grid%toposoil , grid%ht , &
287 ids , ide , jds , jde , 1 , 1 , &
288 ims , ime , jms , jme , 1 , 1 , &
289 ips , ipe , jps , jpe , 1 , 1 )
291 END IF
293 ! Filter the input topography if this is a polar projection.
295 IF ( config_flags%map_proj .EQ. PROJ_CASSINI ) THEN
296 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
298 ! We stick the topo and map fac in an unused 3d array. The map scale
299 ! factor and computational latitude are passed along for the ride
300 ! (part of the transpose process - we only do 3d arrays) to determine
301 ! "how many" values are used to compute the mean. We want a number
302 ! that is consistent with the original grid resolution.
305 DO j = jts, MIN(jte,jde-1)
306 DO k = kts, kte
307 DO i = its, MIN(ite,ide-1)
308 grid%t_init(i,k,j) = 1.
309 END DO
310 END DO
311 DO i = its, MIN(ite,ide-1)
312 grid%t_init(i,1,j) = grid%ht(i,j)
313 grid%t_init(i,2,j) = grid%msftx(i,j)
314 grid%t_init(i,3,j) = grid%clat(i,j)
315 END DO
316 END DO
318 # include "XPOSE_POLAR_FILTER_TOPO_z2x.inc"
320 ! Retrieve the 2d arrays for topo, map factors, and the
321 ! computational latitude.
323 DO j = jpsx, MIN(jpex,jde-1)
324 DO i = ipsx, MIN(ipex,ide-1)
325 grid%ht_xxx(i,j) = grid%t_xxx(i,1,j)
326 grid%mf_xxx(i,j) = grid%t_xxx(i,2,j)
327 grid%clat_xxx(i,j) = grid%t_xxx(i,3,j)
328 END DO
329 END DO
331 ! Get a mean topo field that is consistent with the grid
332 ! distance on each computational latitude loop.
334 CALL filter_topo ( grid%ht_xxx , grid%clat_xxx , grid%mf_xxx , &
335 grid%fft_filter_lat , &
336 ids, ide, jds, jde, 1 , 1 , &
337 imsx, imex, jmsx, jmex, 1, 1, &
338 ipsx, ipex, jpsx, jpex, 1, 1 )
340 ! Stick the filtered topo back into the dummy 3d array to
341 ! transpose it back to "all z on a patch".
343 DO j = jpsx, MIN(jpex,jde-1)
344 DO i = ipsx, MIN(ipex,ide-1)
345 grid%t_xxx(i,1,j) = grid%ht_xxx(i,j)
346 END DO
347 END DO
349 # include "XPOSE_POLAR_FILTER_TOPO_x2z.inc"
351 ! Get the un-transposed topo data.
353 DO j = jts, MIN(jte,jde-1)
354 DO i = its, MIN(ite,ide-1)
355 grid%ht(i,j) = grid%t_init(i,1,j)
356 END DO
357 END DO
358 #else
359 CALL filter_topo ( grid%ht , grid%clat , grid%msftx , &
360 grid%fft_filter_lat , &
361 ids, ide, jds, jde, 1,1, &
362 ims, ime, jms, jme, 1,1, &
363 its, ite, jts, jte, 1,1 )
364 #endif
365 END IF
367 ! If we have any input low-res surface pressure, we store it.
369 IF ( flag_psfc .EQ. 1 ) THEN
370 DO j = jts, MIN(jte,jde-1)
371 DO i = its, MIN(ite,ide-1)
372 grid%psfc_gc(i,j) = grid%psfc(i,j)
373 grid%p_gc(i,1,j) = grid%psfc(i,j)
374 END DO
375 END DO
376 END IF
378 ! If we have the low-resolution surface elevation, stick that in the
379 ! "input" locations of the 3d height. We still have the "hi-res" topo
380 ! stuck in the grid%ht array. The grid%landmask if test is required as some sources
381 ! have ZERO elevation over water (thank you very much).
383 IF ( flag_soilhgt .EQ. 1) THEN
384 DO j = jts, MIN(jte,jde-1)
385 DO i = its, MIN(ite,ide-1)
386 ! IF ( grid%landmask(i,j) .GT. 0.5 ) THEN
387 grid%ght_gc(i,1,j) = grid%toposoil(i,j)
388 grid%ht_gc(i,j)= grid%toposoil(i,j)
389 ! END IF
390 END DO
391 END DO
392 END IF
394 ! Assign surface fields with original input values. If this is hybrid data,
395 ! the values are not exactly representative. However - this is only for
396 ! plotting purposes and such at the 0h of the forecast, so we are not all that
397 ! worried.
399 DO j = jts, min(jde-1,jte)
400 DO i = its, min(ide,ite)
401 grid%u10(i,j)=grid%u_gc(i,1,j)
402 END DO
403 END DO
405 DO j = jts, min(jde,jte)
406 DO i = its, min(ide-1,ite)
407 grid%v10(i,j)=grid%v_gc(i,1,j)
408 END DO
409 END DO
411 DO j = jts, min(jde-1,jte)
412 DO i = its, min(ide-1,ite)
413 grid%t2(i,j)=grid%t_gc(i,1,j)
414 END DO
415 END DO
417 IF ( flag_qv .EQ. 1 ) THEN
418 DO j = jts, min(jde-1,jte)
419 DO i = its, min(ide-1,ite)
420 grid%q2(i,j)=grid%qv_gc(i,1,j)
421 END DO
422 END DO
423 END IF
425 ! The number of vertical levels in the input data. There is no staggering for
426 ! different variables.
428 num_metgrid_levels = grid%num_metgrid_levels
430 ! The requested ptop for real data cases.
432 p_top_requested = grid%p_top_requested
434 ! Compute the top pressure, grid%p_top. For isobaric data, this is just the
435 ! top level. For the generalized vertical coordinate data, we find the
436 ! max pressure on the top level. We have to be careful of two things:
437 ! 1) the value has to be communicated, 2) the value can not increase
438 ! at subsequent times from the initial value.
440 IF ( internal_time_loop .EQ. 1 ) THEN
441 CALL find_p_top ( grid%p_gc , grid%p_top , &
442 ids , ide , jds , jde , 1 , num_metgrid_levels , &
443 ims , ime , jms , jme , 1 , num_metgrid_levels , &
444 its , ite , jts , jte , 1 , num_metgrid_levels )
446 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
447 grid%p_top = wrf_dm_max_real ( grid%p_top )
448 #endif
450 ! Compare the requested grid%p_top with the value available from the input data.
452 IF ( p_top_requested .LT. grid%p_top ) THEN
453 print *,'p_top_requested = ',p_top_requested
454 print *,'allowable grid%p_top in data = ',grid%p_top
455 CALL wrf_error_fatal ( 'p_top_requested < grid%p_top possible from data' )
456 END IF
458 ! The grid%p_top valus is the max of what is available from the data and the
459 ! requested value. We have already compared <, so grid%p_top is directly set to
460 ! the value in the namelist.
462 grid%p_top = p_top_requested
464 ! For subsequent times, we have to remember what the grid%p_top for the first
465 ! time was. Why? If we have a generalized vert coordinate, the grid%p_top value
466 ! could fluctuate.
468 p_top_save = grid%p_top
470 ELSE
471 CALL find_p_top ( grid%p_gc , grid%p_top , &
472 ids , ide , jds , jde , 1 , num_metgrid_levels , &
473 ims , ime , jms , jme , 1 , num_metgrid_levels , &
474 its , ite , jts , jte , 1 , num_metgrid_levels )
476 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
477 grid%p_top = wrf_dm_max_real ( grid%p_top )
478 #endif
479 IF ( grid%p_top .GT. p_top_save ) THEN
480 print *,'grid%p_top from last time period = ',p_top_save
481 print *,'grid%p_top from this time period = ',grid%p_top
482 CALL wrf_error_fatal ( 'grid%p_top > previous value' )
483 END IF
484 grid%p_top = p_top_save
485 ENDIF
487 ! Get the monthly values interpolated to the current date for the traditional monthly
488 ! fields of green-ness fraction and background albedo.
490 CALL monthly_interp_to_date ( grid%greenfrac , current_date , grid%vegfra , &
491 ids , ide , jds , jde , kds , kde , &
492 ims , ime , jms , jme , kms , kme , &
493 its , ite , jts , jte , kts , kte )
495 CALL monthly_interp_to_date ( grid%albedo12m , current_date , grid%albbck , &
496 ids , ide , jds , jde , kds , kde , &
497 ims , ime , jms , jme , kms , kme , &
498 its , ite , jts , jte , kts , kte )
500 ! Get the min/max of each i,j for the monthly green-ness fraction.
502 CALL monthly_min_max ( grid%greenfrac , grid%shdmin , grid%shdmax , &
503 ids , ide , jds , jde , kds , kde , &
504 ims , ime , jms , jme , kms , kme , &
505 its , ite , jts , jte , kts , kte )
507 ! The model expects the green-ness values in percent, not fraction.
509 DO j = jts, MIN(jte,jde-1)
510 DO i = its, MIN(ite,ide-1)
511 grid%vegfra(i,j) = grid%vegfra(i,j) * 100.
512 grid%shdmax(i,j) = grid%shdmax(i,j) * 100.
513 grid%shdmin(i,j) = grid%shdmin(i,j) * 100.
514 END DO
515 END DO
517 ! The model expects the albedo fields as a fraction, not a percent. Set the
518 ! water values to 8%.
520 DO j = jts, MIN(jte,jde-1)
521 DO i = its, MIN(ite,ide-1)
522 grid%albbck(i,j) = grid%albbck(i,j) / 100.
523 grid%snoalb(i,j) = grid%snoalb(i,j) / 100.
524 IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
525 grid%albbck(i,j) = 0.08
526 grid%snoalb(i,j) = 0.08
527 END IF
528 END DO
529 END DO
531 ! Compute the mixing ratio from the input relative humidity.
533 IF ( flag_qv .NE. 1 ) THEN
534 CALL rh_to_mxrat (grid%rh_gc, grid%t_gc, grid%p_gc, grid%qv_gc , .TRUE. , &
535 ids , ide , jds , jde , 1 , num_metgrid_levels , &
536 ims , ime , jms , jme , 1 , num_metgrid_levels , &
537 its , ite , jts , jte , 1 , num_metgrid_levels )
538 END IF
540 ! Two ways to get the surface pressure. 1) If we have the low-res input surface
541 ! pressure and the low-res topography, then we can do a simple hydrostatic
542 ! relation. 2) Otherwise we compute the surface pressure from the sea-level
543 ! pressure.
544 ! Note that on output, grid%psfc is now hi-res. The low-res surface pressure and
545 ! elevation are grid%psfc_gc and grid%ht_gc (same as grid%ght_gc(k=1)).
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 , tiso , &
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 ! Set the temperature of the inland lakes to tavgsfc if the temperature is available
853 ! and islake is > num_veg_cat
855 num_veg_cat = SIZE ( grid%landusef , DIM=2 )
856 CALL nl_get_iswater ( grid%id , grid%iswater )
857 CALL nl_get_islake ( grid%id , grid%islake )
859 IF ( grid%islake < 0 ) THEN
860 CALL wrf_debug ( 0 , 'Old data, no inland lake information')
861 ELSE
862 IF ( we_have_tavgsfc ) THEN
864 CALL wrf_debug ( 0 , 'Using inland lakes with average surface temperature')
865 DO j=jts,MIN(jde-1,jte)
866 DO i=its,MIN(ide-1,ite)
867 IF ( grid%landusef(i,grid%islake,j) >= 0.5 ) THEN
868 grid%sst(i,j) = grid%tavgsfc(i,j)
869 grid%tsk(i,j) = grid%tavgsfc(i,j)
870 END IF
871 END DO
872 END DO
874 ELSE ! We don't have tavgsfc
876 CALL wrf_debug ( 0 , 'No average surface temperature for use with inland lakes')
878 END IF
879 DO j=jts,MIN(jde-1,jte)
880 DO i=its,MIN(ide-1,ite)
881 grid%landusef(i,grid%iswater,j) = grid%landusef(i,grid%iswater,j) + &
882 grid%landusef(i,grid%islake,j)
883 grid%landusef(i,grid%islake,j) = 0.
884 END DO
885 END DO
887 END IF
889 ! Save the grid%tsk field for later use in the sea ice surface temperature
890 ! for the Noah LSM scheme.
892 DO j = jts, MIN(jte,jde-1)
893 DO i = its, MIN(ite,ide-1)
894 grid%tsk_save(i,j) = grid%tsk(i,j)
895 END DO
896 END DO
898 ! Protect against bad grid%tsk values over water by supplying grid%sst (if it is
899 ! available, and if the grid%sst is reasonable).
901 DO j = jts, MIN(jde-1,jte)
902 DO i = its, MIN(ide-1,ite)
903 IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
904 ( grid%sst(i,j) .GT. 170. ) .AND. ( grid%sst(i,j) .LT. 400. ) ) THEN
905 grid%tsk(i,j) = grid%sst(i,j)
906 ENDIF
907 END DO
908 END DO
910 ! Take the data from the input file and store it in the variables that
911 ! use the WRF naming and ordering conventions.
913 DO j = jts, MIN(jte,jde-1)
914 DO i = its, MIN(ite,ide-1)
915 IF ( grid%snow(i,j) .GE. 10. ) then
916 grid%snowc(i,j) = 1.
917 ELSE
918 grid%snowc(i,j) = 0.0
919 END IF
920 END DO
921 END DO
923 ! Set flag integers for presence of snowh and soilw fields
925 grid%ifndsnowh = flag_snowh
926 IF (num_sw_levels_input .GE. 1) THEN
927 grid%ifndsoilw = 1
928 ELSE
929 grid%ifndsoilw = 0
930 END IF
932 ! We require input data for the various LSM schemes.
934 enough_data : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
936 CASE (LSMSCHEME)
937 IF ( num_st_levels_input .LT. 2 ) THEN
938 CALL wrf_error_fatal ( 'Not enough soil temperature data for Noah LSM scheme.')
939 END IF
941 CASE (RUCLSMSCHEME)
942 IF ( num_st_levels_input .LT. 2 ) THEN
943 CALL wrf_error_fatal ( 'Not enough soil temperature data for RUC LSM scheme.')
944 END IF
946 CASE (PXLSMSCHEME)
947 IF ( num_st_levels_input .LT. 2 ) THEN
948 CALL wrf_error_fatal ( 'Not enough soil temperature data for P-X LSM scheme.')
949 END IF
951 END SELECT enough_data
953 interpolate_soil_tmw : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
955 CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME, PXLSMSCHEME )
956 CALL process_soil_real ( grid%tsk , grid%tmn , grid%tavgsfc, &
957 grid%landmask , grid%sst , grid%ht, grid%toposoil, &
958 st_input , sm_input , sw_input , &
959 st_levels_input , sm_levels_input , sw_levels_input , &
960 grid%zs , grid%dzs , grid%tslb , grid%smois , grid%sh2o , &
961 flag_sst , flag_tavgsfc, &
962 flag_soilhgt, flag_soil_layers, flag_soil_levels, &
963 ids , ide , jds , jde , kds , kde , &
964 ims , ime , jms , jme , kms , kme , &
965 its , ite , jts , jte , kts , kte , &
966 model_config_rec%sf_surface_physics(grid%id) , &
967 model_config_rec%num_soil_layers , &
968 model_config_rec%real_data_init_type , &
969 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
970 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
972 END SELECT interpolate_soil_tmw
974 ! Adjustments for the seaice field PRIOR to the grid%tslb computations. This is
975 ! is for the 5-layer scheme.
977 num_veg_cat = SIZE ( grid%landusef , DIM=2 )
978 num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
979 num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
980 CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
981 CALL nl_get_isice ( grid%id , grid%isice )
982 CALL nl_get_iswater ( grid%id , grid%iswater )
983 CALL adjust_for_seaice_pre ( grid%xice , grid%landmask , grid%tsk , grid%ivgtyp , grid%vegcat , grid%lu_index , &
984 grid%xland , grid%landusef , grid%isltyp , grid%soilcat , grid%soilctop , &
985 grid%soilcbot , grid%tmn , &
986 grid%seaice_threshold , &
987 config_flags%fractional_seaice, &
988 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
989 grid%iswater , grid%isice , &
990 model_config_rec%sf_surface_physics(grid%id) , &
991 ids , ide , jds , jde , kds , kde , &
992 ims , ime , jms , jme , kms , kme , &
993 its , ite , jts , jte , kts , kte )
995 ! surface_input_source=1 => use data from static file (fractional category as input)
996 ! surface_input_source=2 => use data from grib file (dominant category as input)
998 IF ( config_flags%surface_input_source .EQ. 1 ) THEN
999 grid%vegcat (its,jts) = 0
1000 grid%soilcat(its,jts) = 0
1001 END IF
1003 ! Generate the vegetation and soil category information from the fractional input
1004 ! data, or use the existing dominant category fields if they exist.
1006 IF ( ( grid%soilcat(its,jts) .LT. 0.5 ) .AND. ( grid%vegcat(its,jts) .LT. 0.5 ) ) THEN
1008 num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1009 num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1010 num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1012 CALL process_percent_cat_new ( grid%landmask , &
1013 grid%landusef , grid%soilctop , grid%soilcbot , &
1014 grid%isltyp , grid%ivgtyp , &
1015 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1016 ids , ide , jds , jde , kds , kde , &
1017 ims , ime , jms , jme , kms , kme , &
1018 its , ite , jts , jte , kts , kte , &
1019 model_config_rec%iswater(grid%id) )
1021 ! Make all the veg/soil parms the same so as not to confuse the developer.
1023 DO j = jts , MIN(jde-1,jte)
1024 DO i = its , MIN(ide-1,ite)
1025 grid%vegcat(i,j) = grid%ivgtyp(i,j)
1026 grid%soilcat(i,j) = grid%isltyp(i,j)
1027 END DO
1028 END DO
1030 ELSE
1032 ! Do we have dominant soil and veg data from the input already?
1034 IF ( grid%soilcat(its,jts) .GT. 0.5 ) THEN
1035 DO j = jts, MIN(jde-1,jte)
1036 DO i = its, MIN(ide-1,ite)
1037 grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
1038 END DO
1039 END DO
1040 END IF
1041 IF ( grid%vegcat(its,jts) .GT. 0.5 ) THEN
1042 DO j = jts, MIN(jde-1,jte)
1043 DO i = its, MIN(ide-1,ite)
1044 grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
1045 END DO
1046 END DO
1047 END IF
1049 END IF
1051 ! Land use assignment.
1053 DO j = jts, MIN(jde-1,jte)
1054 DO i = its, MIN(ide-1,ite)
1055 grid%lu_index(i,j) = grid%ivgtyp(i,j)
1056 IF ( grid%lu_index(i,j) .NE. model_config_rec%iswater(grid%id) ) THEN
1057 grid%landmask(i,j) = 1
1058 grid%xland(i,j) = 1
1059 ELSE
1060 grid%landmask(i,j) = 0
1061 grid%xland(i,j) = 2
1062 END IF
1063 END DO
1064 END DO
1066 ! Fix grid%tmn and grid%tsk.
1068 fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1070 CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME, PXLSMSCHEME )
1071 DO j = jts, MIN(jde-1,jte)
1072 DO i = its, MIN(ide-1,ite)
1073 IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
1074 ( grid%sst(i,j) .GT. 170. ) .AND. ( grid%sst(i,j) .LT. 400. ) ) THEN
1075 grid%tmn(i,j) = grid%sst(i,j)
1076 grid%tsk(i,j) = grid%sst(i,j)
1077 ELSE IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
1078 grid%tmn(i,j) = grid%tsk(i,j)
1079 END IF
1080 END DO
1081 END DO
1082 END SELECT fix_tsk_tmn
1084 ! Is the grid%tsk reasonable?
1086 IF ( internal_time_loop .NE. 1 ) THEN
1087 DO j = jts, MIN(jde-1,jte)
1088 DO i = its, MIN(ide-1,ite)
1089 IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
1090 grid%tsk(i,j) = grid%t_2(i,1,j)
1091 END IF
1092 END DO
1093 END DO
1094 ELSE
1095 DO j = jts, MIN(jde-1,jte)
1096 DO i = its, MIN(ide-1,ite)
1097 IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
1098 print *,'error in the grid%tsk'
1099 print *,'i,j=',i,j
1100 print *,'grid%landmask=',grid%landmask(i,j)
1101 print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1102 if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
1103 grid%tsk(i,j)=grid%tmn(i,j)
1104 else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1105 grid%tsk(i,j)=grid%sst(i,j)
1106 else
1107 CALL wrf_error_fatal ( 'grid%tsk unreasonable' )
1108 end if
1109 END IF
1110 END DO
1111 END DO
1112 END IF
1114 ! Is the grid%tmn reasonable?
1116 DO j = jts, MIN(jde-1,jte)
1117 DO i = its, MIN(ide-1,ite)
1118 IF ( ( ( grid%tmn(i,j) .LT. 170. ) .OR. ( grid%tmn(i,j) .GT. 400. ) ) &
1119 .AND. ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
1120 IF ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) THEN
1121 print *,'error in the grid%tmn'
1122 print *,'i,j=',i,j
1123 print *,'grid%landmask=',grid%landmask(i,j)
1124 print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1125 END IF
1127 if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
1128 grid%tmn(i,j)=grid%tsk(i,j)
1129 else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1130 grid%tmn(i,j)=grid%sst(i,j)
1131 else
1132 CALL wrf_error_fatal ( 'grid%tmn unreasonable' )
1133 endif
1134 END IF
1135 END DO
1136 END DO
1139 ! Minimum soil values, residual, from RUC LSM scheme. For input from Noah or EC, and using
1140 ! RUC LSM scheme, this must be subtracted from the input total soil moisture. For
1141 ! input RUC data and using the Noah LSM scheme, this value must be added to the soil
1142 ! moisture input.
1144 lqmi(1:num_soil_top_cat) = &
1145 (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
1146 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
1147 0.004, 0.065 /)
1148 ! 0.004, 0.065, 0.020, 0.004, 0.008 /) ! has extra levels for playa, lava, and white sand
1150 ! At the initial time we care about values of soil moisture and temperature, other times are
1151 ! ignored by the model, so we ignore them, too.
1153 IF ( domain_ClockIsStartTime(grid) ) THEN
1154 account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1156 CASE ( LSMSCHEME )
1157 iicount = 0
1158 IF ( flag_soil_layers == 1 ) THEN
1159 DO j = jts, MIN(jde-1,jte)
1160 DO i = its, MIN(ide-1,ite)
1161 IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 170 ) .and. &
1162 ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
1163 print *,'Noah -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
1164 iicount = iicount + 1
1165 grid%smois(i,:,j) = 0.005
1166 END IF
1167 END DO
1168 END DO
1169 IF ( iicount .GT. 0 ) THEN
1170 print *,'Noah -> Noah: total number of small soil moisture locations = ',iicount
1171 END IF
1172 ELSE IF ( flag_soil_levels == 1 ) THEN
1173 DO j = jts, MIN(jde-1,jte)
1174 DO i = its, MIN(ide-1,ite)
1175 grid%smois(i,:,j) = grid%smois(i,:,j) + lqmi(grid%isltyp(i,j))
1176 END DO
1177 END DO
1178 DO j = jts, MIN(jde-1,jte)
1179 DO i = its, MIN(ide-1,ite)
1180 IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 170 ) .and. &
1181 ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
1182 print *,'RUC -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
1183 iicount = iicount + 1
1184 grid%smois(i,:,j) = 0.005
1185 END IF
1186 END DO
1187 END DO
1188 IF ( iicount .GT. 0 ) THEN
1189 print *,'RUC -> Noah: total number of small soil moisture locations = ',iicount
1190 END IF
1191 END IF
1193 CASE ( RUCLSMSCHEME )
1194 iicount = 0
1195 IF ( flag_soil_layers == 1 ) THEN
1196 DO j = jts, MIN(jde-1,jte)
1197 DO i = its, MIN(ide-1,ite)
1198 grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) - lqmi(grid%isltyp(i,j)) , 0. )
1199 END DO
1200 END DO
1201 ELSE IF ( flag_soil_levels == 1 ) THEN
1202 ! no op
1203 END IF
1205 CASE ( PXLSMSCHEME )
1206 iicount = 0
1207 IF ( flag_soil_layers == 1 ) THEN
1208 DO j = jts, MIN(jde-1,jte)
1209 DO i = its, MIN(ide-1,ite)
1210 grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) - lqmi(grid%isltyp(i,j)) , 0. )
1211 END DO
1212 END DO
1213 ELSE IF ( flag_soil_levels == 1 ) THEN
1214 ! no op
1215 END IF
1217 END SELECT account_for_zero_soil_moisture
1218 END IF
1220 ! Is the grid%tslb reasonable?
1222 IF ( internal_time_loop .NE. 1 ) THEN
1223 DO j = jts, MIN(jde-1,jte)
1224 DO ns = 1 , model_config_rec%num_soil_layers
1225 DO i = its, MIN(ide-1,ite)
1226 IF ( grid%tslb(i,ns,j) .LT. 170 .or. grid%tslb(i,ns,j) .GT. 400. ) THEN
1227 grid%tslb(i,ns,j) = grid%t_2(i,1,j)
1228 grid%smois(i,ns,j) = 0.3
1229 END IF
1230 END DO
1231 END DO
1232 END DO
1233 ELSE
1234 DO j = jts, MIN(jde-1,jte)
1235 DO i = its, MIN(ide-1,ite)
1236 IF ( ( ( grid%tslb(i,1,j) .LT. 170. ) .OR. ( grid%tslb(i,1,j) .GT. 400. ) ) .AND. &
1237 ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
1238 IF ( ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) .AND. &
1239 ( model_config_rec%sf_surface_physics(grid%id) .NE. RUCLSMSCHEME ).AND. &
1240 ( model_config_rec%sf_surface_physics(grid%id) .NE. PXLSMSCHEME ) ) THEN
1241 print *,'error in the grid%tslb'
1242 print *,'i,j=',i,j
1243 print *,'grid%landmask=',grid%landmask(i,j)
1244 print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1245 print *,'grid%tslb = ',grid%tslb(i,:,j)
1246 print *,'old grid%smois = ',grid%smois(i,:,j)
1247 grid%smois(i,1,j) = 0.3
1248 grid%smois(i,2,j) = 0.3
1249 grid%smois(i,3,j) = 0.3
1250 grid%smois(i,4,j) = 0.3
1251 END IF
1253 IF ( (grid%tsk(i,j).GT.170. .AND. grid%tsk(i,j).LT.400.) .AND. &
1254 (grid%tmn(i,j).GT.170. .AND. grid%tmn(i,j).LT.400.) ) THEN
1255 fake_soil_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1256 CASE ( SLABSCHEME )
1257 DO ns = 1 , model_config_rec%num_soil_layers
1258 grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
1259 grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
1260 END DO
1261 CASE ( LSMSCHEME , RUCLSMSCHEME, PXLSMSCHEME )
1262 CALL wrf_error_fatal ( 'Assigning constant soil moisture, bad idea')
1263 DO ns = 1 , model_config_rec%num_soil_layers
1264 grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
1265 grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
1266 END DO
1267 END SELECT fake_soil_temp
1268 else if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
1269 CALL wrf_error_fatal ( 'grid%tslb unreasonable 1' )
1270 DO ns = 1 , model_config_rec%num_soil_layers
1271 grid%tslb(i,ns,j)=grid%tsk(i,j)
1272 END DO
1273 else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1274 CALL wrf_error_fatal ( 'grid%tslb unreasonable 2' )
1275 DO ns = 1 , model_config_rec%num_soil_layers
1276 grid%tslb(i,ns,j)=grid%sst(i,j)
1277 END DO
1278 else if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
1279 CALL wrf_error_fatal ( 'grid%tslb unreasonable 3' )
1280 DO ns = 1 , model_config_rec%num_soil_layers
1281 grid%tslb(i,ns,j)=grid%tmn(i,j)
1282 END DO
1283 else
1284 CALL wrf_error_fatal ( 'grid%tslb unreasonable 4' )
1285 endif
1286 END IF
1287 END DO
1288 END DO
1289 END IF
1291 ! Adjustments for the seaice field AFTER the grid%tslb computations. This is
1292 ! is for the Noah LSM scheme.
1294 num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1295 num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1296 num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1297 CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
1298 CALL nl_get_isice ( grid%id , grid%isice )
1299 CALL nl_get_iswater ( grid%id , grid%iswater )
1300 CALL adjust_for_seaice_post ( grid%xice , grid%landmask , grid%tsk , grid%tsk_save , &
1301 grid%ivgtyp , grid%vegcat , grid%lu_index , &
1302 grid%xland , grid%landusef , grid%isltyp , grid%soilcat , &
1303 grid%soilctop , &
1304 grid%soilcbot , grid%tmn , grid%vegfra , &
1305 grid%tslb , grid%smois , grid%sh2o , &
1306 grid%seaice_threshold , &
1307 config_flags%fractional_seaice, &
1308 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1309 model_config_rec%num_soil_layers , &
1310 grid%iswater , grid%isice , &
1311 model_config_rec%sf_surface_physics(grid%id) , &
1312 ids , ide , jds , jde , kds , kde , &
1313 ims , ime , jms , jme , kms , kme , &
1314 its , ite , jts , jte , kts , kte )
1316 ! Let us make sure (again) that the grid%landmask and the veg/soil categories match.
1318 oops1=0
1319 oops2=0
1320 DO j = jts, MIN(jde-1,jte)
1321 DO i = its, MIN(ide-1,ite)
1322 IF ( ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. &
1323 ( grid%ivgtyp(i,j) .NE. config_flags%iswater .OR. grid%isltyp(i,j) .NE. 14 ) ) .OR. &
1324 ( ( grid%landmask(i,j) .GT. 0.5 ) .AND. &
1325 ( grid%ivgtyp(i,j) .EQ. config_flags%iswater .OR. grid%isltyp(i,j) .EQ. 14 ) ) ) THEN
1326 IF ( grid%tslb(i,1,j) .GT. 1. ) THEN
1327 oops1=oops1+1
1328 grid%ivgtyp(i,j) = 5
1329 grid%isltyp(i,j) = 8
1330 grid%landmask(i,j) = 1
1331 grid%xland(i,j) = 1
1332 ELSE IF ( grid%sst(i,j) .GT. 1. ) THEN
1333 oops2=oops2+1
1334 grid%ivgtyp(i,j) = config_flags%iswater
1335 grid%isltyp(i,j) = 14
1336 grid%landmask(i,j) = 0
1337 grid%xland(i,j) = 2
1338 ELSE
1339 print *,'the grid%landmask and soil/veg cats do not match'
1340 print *,'i,j=',i,j
1341 print *,'grid%landmask=',grid%landmask(i,j)
1342 print *,'grid%ivgtyp=',grid%ivgtyp(i,j)
1343 print *,'grid%isltyp=',grid%isltyp(i,j)
1344 print *,'iswater=', config_flags%iswater
1345 print *,'grid%tslb=',grid%tslb(i,:,j)
1346 print *,'grid%sst=',grid%sst(i,j)
1347 CALL wrf_error_fatal ( 'mismatch_landmask_ivgtyp' )
1348 END IF
1349 END IF
1350 END DO
1351 END DO
1352 if (oops1.gt.0) then
1353 print *,'points artificially set to land : ',oops1
1354 endif
1355 if(oops2.gt.0) then
1356 print *,'points artificially set to water: ',oops2
1357 endif
1358 ! fill grid%sst array with grid%tsk if missing in real input (needed for time-varying grid%sst in wrf)
1359 DO j = jts, MIN(jde-1,jte)
1360 DO i = its, MIN(ide-1,ite)
1361 IF ( flag_sst .NE. 1 ) THEN
1362 grid%sst(i,j) = grid%tsk(i,j)
1363 ENDIF
1364 END DO
1365 END DO
1367 ! From the full level data, we can get the half levels, reciprocals, and layer
1368 ! thicknesses. These are all defined at half level locations, so one less level.
1369 ! We allow the vertical coordinate to *accidently* come in upside down. We want
1370 ! the first full level to be the ground surface.
1372 ! Check whether grid%znw (full level) data are truly full levels. If not, we need to adjust them
1373 ! to be full levels.
1374 ! in this test, we check if grid%znw(1) is neither 0 nor 1 (within a tolerance of 10**-5)
1376 were_bad = .false.
1377 IF ( ( (grid%znw(1).LT.(1-1.E-5) ) .OR. ( grid%znw(1).GT.(1+1.E-5) ) ).AND. &
1378 ( (grid%znw(1).LT.(0-1.E-5) ) .OR. ( grid%znw(1).GT.(0+1.E-5) ) ) ) THEN
1379 were_bad = .true.
1380 print *,'Your grid%znw input values are probably half-levels. '
1381 print *,grid%znw
1382 print *,'WRF expects grid%znw values to be full levels. '
1383 print *,'Adjusting now to full levels...'
1384 ! We want to ignore the first value if it's negative
1385 IF (grid%znw(1).LT.0) THEN
1386 grid%znw(1)=0
1387 END IF
1388 DO k=2,kde
1389 grid%znw(k)=2*grid%znw(k)-grid%znw(k-1)
1390 END DO
1391 END IF
1393 ! Let's check our changes
1395 IF ( ( ( grid%znw(1) .LT. (1-1.E-5) ) .OR. ( grid%znw(1) .GT. (1+1.E-5) ) ).AND. &
1396 ( ( grid%znw(1) .LT. (0-1.E-5) ) .OR. ( grid%znw(1) .GT. (0+1.E-5) ) ) ) THEN
1397 print *,'The input grid%znw height values were half-levels or erroneous. '
1398 print *,'Attempts to treat the values as half-levels and change them '
1399 print *,'to valid full levels failed.'
1400 CALL wrf_error_fatal("bad grid%znw values from input files")
1401 ELSE IF ( were_bad ) THEN
1402 print *,'...adjusted. grid%znw array now contains full eta level values. '
1403 ENDIF
1405 IF ( grid%znw(1) .LT. grid%znw(kde) ) THEN
1406 DO k=1, kde/2
1407 hold_znw = grid%znw(k)
1408 grid%znw(k)=grid%znw(kde+1-k)
1409 grid%znw(kde+1-k)=hold_znw
1410 END DO
1411 END IF
1413 DO k=1, kde-1
1414 grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
1415 grid%rdnw(k) = 1./grid%dnw(k)
1416 grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
1417 END DO
1419 ! Now the same sort of computations with the half eta levels, even ANOTHER
1420 ! level less than the one above.
1422 DO k=2, kde-1
1423 grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
1424 grid%rdn(k) = 1./grid%dn(k)
1425 grid%fnp(k) = .5* grid%dnw(k )/grid%dn(k)
1426 grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
1427 END DO
1429 ! Scads of vertical coefficients.
1431 cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2)
1432 cof2 = grid%dn(2) /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
1434 grid%cf1 = grid%fnp(2) + cof1
1435 grid%cf2 = grid%fnm(2) - cof1 - cof2
1436 grid%cf3 = cof2
1438 grid%cfn = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
1439 grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
1441 ! Inverse grid distances.
1443 grid%rdx = 1./config_flags%dx
1444 grid%rdy = 1./config_flags%dy
1446 ! Some of the many weird geopotential initializations that we'll see today: grid%ph0 is total,
1447 ! and grid%ph_2 is a perturbation from the base state geopotential. We set the base geopotential
1448 ! at the lowest level to terrain elevation * gravity.
1450 DO j=jts,jte
1451 DO i=its,ite
1452 grid%ph0(i,1,j) = grid%ht(i,j) * g
1453 grid%ph_2(i,1,j) = 0.
1454 END DO
1455 END DO
1457 ! Base state potential temperature and inverse density (alpha = 1/rho) from
1458 ! the half eta levels and the base-profile surface pressure. Compute 1/rho
1459 ! from equation of state. The potential temperature is a perturbation from t0.
1461 DO j = jts, MIN(jte,jde-1)
1462 DO i = its, MIN(ite,ide-1)
1464 ! Base state pressure is a function of eta level and terrain, only, plus
1465 ! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
1466 ! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
1468 p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
1471 DO k = 1, kte-1
1472 grid%php(i,k,j) = grid%znw(k)*(p_surf - grid%p_top) + grid%p_top ! temporary, full lev base pressure
1473 grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
1474 temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
1475 ! temp = t00 + A*LOG(grid%pb(i,k,j)/p00)
1476 grid%t_init(i,k,j) = temp*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
1477 grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
1478 END DO
1480 ! Base state mu is defined as base state surface pressure minus grid%p_top
1482 grid%mub(i,j) = p_surf - grid%p_top
1484 ! Dry surface pressure is defined as the following (this mu is from the input file
1485 ! computed from the dry pressure). Here the dry pressure is just reconstituted.
1487 pd_surf = grid%mu0(i,j) + grid%p_top
1489 ! Integrate base geopotential, starting at terrain elevation. This assures that
1490 ! the base state is in exact hydrostatic balance with respect to the model equations.
1491 ! This field is on full levels.
1493 grid%phb(i,1,j) = grid%ht(i,j) * g
1494 DO k = 2,kte
1495 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)
1496 END DO
1497 END DO
1498 END DO
1500 ! Fill in the outer rows and columns to allow us to be sloppy.
1502 IF ( ite .EQ. ide ) THEN
1503 i = ide
1504 DO j = jts, MIN(jde-1,jte)
1505 grid%mub(i,j) = grid%mub(i-1,j)
1506 grid%mu_2(i,j) = grid%mu_2(i-1,j)
1507 DO k = 1, kte-1
1508 grid%pb(i,k,j) = grid%pb(i-1,k,j)
1509 grid%t_init(i,k,j) = grid%t_init(i-1,k,j)
1510 grid%alb(i,k,j) = grid%alb(i-1,k,j)
1511 END DO
1512 DO k = 1, kte
1513 grid%phb(i,k,j) = grid%phb(i-1,k,j)
1514 END DO
1515 END DO
1516 END IF
1518 IF ( jte .EQ. jde ) THEN
1519 j = jde
1520 DO i = its, ite
1521 grid%mub(i,j) = grid%mub(i,j-1)
1522 grid%mu_2(i,j) = grid%mu_2(i,j-1)
1523 DO k = 1, kte-1
1524 grid%pb(i,k,j) = grid%pb(i,k,j-1)
1525 grid%t_init(i,k,j) = grid%t_init(i,k,j-1)
1526 grid%alb(i,k,j) = grid%alb(i,k,j-1)
1527 END DO
1528 DO k = 1, kte
1529 grid%phb(i,k,j) = grid%phb(i,k,j-1)
1530 END DO
1531 END DO
1532 END IF
1534 ! Compute the perturbation dry pressure (grid%mub + grid%mu_2 + ptop = dry grid%psfc).
1536 DO j = jts, min(jde-1,jte)
1537 DO i = its, min(ide-1,ite)
1538 grid%mu_2(i,j) = grid%mu0(i,j) - grid%mub(i,j)
1539 END DO
1540 END DO
1542 ! Fill in the outer rows and columns to allow us to be sloppy.
1544 IF ( ite .EQ. ide ) THEN
1545 i = ide
1546 DO j = jts, MIN(jde-1,jte)
1547 grid%mu_2(i,j) = grid%mu_2(i-1,j)
1548 END DO
1549 END IF
1551 IF ( jte .EQ. jde ) THEN
1552 j = jde
1553 DO i = its, ite
1554 grid%mu_2(i,j) = grid%mu_2(i,j-1)
1555 END DO
1556 END IF
1558 lev500 = 0
1559 DO j = jts, min(jde-1,jte)
1560 DO i = its, min(ide-1,ite)
1562 ! Assign the potential temperature (perturbation from t0) and qv on all the mass
1563 ! point locations.
1565 DO k = 1 , kde-1
1566 grid%t_2(i,k,j) = grid%t_2(i,k,j) - t0
1567 END DO
1569 dpmu = 10001.
1570 loop_count = 0
1572 DO WHILE ( ( ABS(dpmu) .GT. 10. ) .AND. &
1573 ( loop_count .LT. 5 ) )
1575 loop_count = loop_count + 1
1577 ! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
1578 ! equation) down from the top to get the pressure perturbation. First get the pressure
1579 ! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
1581 k = kte-1
1583 qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
1584 qvf2 = 1./(1.+qvf1)
1585 qvf1 = qvf1*qvf2
1587 grid%p(i,k,j) = - 0.5*(grid%mu_2(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
1588 qvf = 1. + rvovrd*moist(i,k,j,P_QV)
1589 grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf&
1590 *(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
1591 grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
1593 ! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
1594 ! inverse density fields (total and perturbation).
1596 DO k=kte-2,1,-1
1597 qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
1598 qvf2 = 1./(1.+qvf1)
1599 qvf1 = qvf1*qvf2
1600 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)
1601 qvf = 1. + rvovrd*moist(i,k,j,P_QV)
1602 grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf* &
1603 (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
1604 grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
1605 END DO
1607 #if 1
1608 ! This is the hydrostatic equation used in the model after the small timesteps. In
1609 ! the model, grid%al (inverse density) is computed from the geopotential.
1611 DO k = 2,kte
1612 grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
1613 grid%dnw(k-1) * ( (grid%mub(i,j)+grid%mu_2(i,j))*grid%al(i,k-1,j) &
1614 + grid%mu_2(i,j)*grid%alb(i,k-1,j) )
1615 grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
1616 END DO
1617 #else
1618 ! Get the perturbation geopotential from the 3d height array from WPS.
1620 DO k = 2,kte
1621 grid%ph_2(i,k,j) = grid%ph0(i,k,j)*g - grid%phb(i,k,j)
1622 END DO
1623 #endif
1625 ! Adjust the column pressure so that the computed 500 mb height is close to the
1626 ! input value (of course, not when we are doing hybrid input).
1628 IF ( ( flag_metgrid .EQ. 1 ) .AND. ( i .EQ. its ) .AND. ( j .EQ. jts ) ) THEN
1629 DO k = 1 , num_metgrid_levels
1630 IF ( ABS ( grid%p_gc(i,k,j) - 50000. ) .LT. 1. ) THEN
1631 lev500 = k
1632 EXIT
1633 END IF
1634 END DO
1635 END IF
1637 ! We only do the adjustment of height if we have the input data on pressure
1638 ! surfaces, and folks have asked to do this option.
1640 IF ( ( flag_metgrid .EQ. 1 ) .AND. &
1641 ( config_flags%adjust_heights ) .AND. &
1642 ( lev500 .NE. 0 ) ) THEN
1644 DO k = 2 , kte-1
1646 ! Get the pressures on the full eta levels (grid%php is defined above as
1647 ! the full-lev base pressure, an easy array to use for 3d space).
1649 pl = grid%php(i,k ,j) + &
1650 ( grid%p(i,k-1 ,j) * ( grid%znw(k ) - grid%znu(k ) ) + &
1651 grid%p(i,k ,j) * ( grid%znu(k-1 ) - grid%znw(k ) ) ) / &
1652 ( grid%znu(k-1 ) - grid%znu(k ) )
1653 pu = grid%php(i,k+1,j) + &
1654 ( grid%p(i,k-1+1,j) * ( grid%znw(k +1) - grid%znu(k+1) ) + &
1655 grid%p(i,k +1,j) * ( grid%znu(k-1+1) - grid%znw(k+1) ) ) / &
1656 ( grid%znu(k-1+1) - grid%znu(k+1) )
1658 ! If these pressure levels trap 500 mb, use them to interpolate
1659 ! to the 500 mb level of the computed height.
1661 IF ( ( pl .GE. 50000. ) .AND. ( pu .LT. 50000. ) ) THEN
1662 zl = ( grid%ph_2(i,k ,j) + grid%phb(i,k ,j) ) / g
1663 zu = ( grid%ph_2(i,k+1,j) + grid%phb(i,k+1,j) ) / g
1665 z500 = ( zl * ( LOG(50000.) - LOG(pu ) ) + &
1666 zu * ( LOG(pl ) - LOG(50000.) ) ) / &
1667 ( LOG(pl) - LOG(pu) )
1668 ! z500 = ( zl * ( (50000.) - (pu ) ) + &
1669 ! zu * ( (pl ) - (50000.) ) ) / &
1670 ! ( (pl) - (pu) )
1672 ! Compute the difference of the 500 mb heights (computed minus input), and
1673 ! then the change in grid%mu_2. The grid%php is still full-levels, base pressure.
1675 dz500 = z500 - grid%ght_gc(i,lev500,j)
1676 tvsfc = ((grid%t_2(i,1,j)+t0)*((grid%p(i,1,j)+grid%php(i,1,j))/p1000mb)**(r_d/cp)) * &
1677 (1.+0.6*moist(i,1,j,P_QV))
1678 dpmu = ( grid%php(i,1,j) + grid%p(i,1,j) ) * EXP ( g * dz500 / ( r_d * tvsfc ) )
1679 dpmu = dpmu - ( grid%php(i,1,j) + grid%p(i,1,j) )
1680 grid%mu_2(i,j) = grid%mu_2(i,j) - dpmu
1681 EXIT
1682 END IF
1684 END DO
1685 ELSE
1686 dpmu = 0.
1687 END IF
1689 END DO
1691 END DO
1692 END DO
1694 ! If this is data from the SI, then we probably do not have the original
1695 ! surface data laying around. Note that these are all the lowest levels
1696 ! of the respective 3d arrays. For surface pressure, we assume that the
1697 ! vertical gradient of grid%p prime is zilch. This is not all that important.
1698 ! These are filled in so that the various plotting routines have something
1699 ! to play with at the initial time for the model.
1701 IF ( flag_metgrid .NE. 1 ) THEN
1702 DO j = jts, min(jde-1,jte)
1703 DO i = its, min(ide,ite)
1704 grid%u10(i,j)=grid%u_2(i,1,j)
1705 END DO
1706 END DO
1708 DO j = jts, min(jde,jte)
1709 DO i = its, min(ide-1,ite)
1710 grid%v10(i,j)=grid%v_2(i,1,j)
1711 END DO
1712 END DO
1714 DO j = jts, min(jde-1,jte)
1715 DO i = its, min(ide-1,ite)
1716 p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
1717 grid%psfc(i,j)=p_surf + grid%p(i,1,j)
1718 grid%q2(i,j)=moist(i,1,j,P_QV)
1719 grid%th2(i,j)=grid%t_2(i,1,j)+300.
1720 grid%t2(i,j)=grid%th2(i,j)*(((grid%p(i,1,j)+grid%pb(i,1,j))/p00)**(r_d/cp))
1721 END DO
1722 END DO
1724 ! If this data is from WPS, then we have previously assigned the surface
1725 ! data for u, v, and t. If we have an input qv, welp, we assigned that one,
1726 ! too. Now we pick up the left overs, and if RH came in - we assign the
1727 ! mixing ratio.
1729 ELSE IF ( flag_metgrid .EQ. 1 ) THEN
1731 DO j = jts, min(jde-1,jte)
1732 DO i = its, min(ide-1,ite)
1733 p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
1734 grid%psfc(i,j)=p_surf + grid%p(i,1,j)
1735 grid%th2(i,j)=grid%t2(i,j)*(p00/(grid%p(i,1,j)+grid%pb(i,1,j)))**(r_d/cp)
1736 END DO
1737 END DO
1738 IF ( flag_qv .NE. 1 ) THEN
1739 DO j = jts, min(jde-1,jte)
1740 DO i = its, min(ide-1,ite)
1741 grid%q2(i,j)=moist(i,1,j,P_QV)
1742 END DO
1743 END DO
1744 END IF
1746 END IF
1748 ! Set flag to denote that we are saving original values of HT, MUB, and
1749 ! PHB for 2-way nesting and cycling.
1751 grid%save_topo_from_real=1
1753 ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
1754 #ifdef DM_PARALLEL
1755 # include "HALO_EM_INIT_1.inc"
1756 # include "HALO_EM_INIT_2.inc"
1757 # include "HALO_EM_INIT_3.inc"
1758 # include "HALO_EM_INIT_4.inc"
1759 # include "HALO_EM_INIT_5.inc"
1760 #endif
1762 RETURN
1764 END SUBROUTINE init_domain_rk
1766 !---------------------------------------------------------------------
1768 SUBROUTINE const_module_initialize ( p00 , t00 , a , tiso )
1769 USE module_configure
1770 IMPLICIT NONE
1771 ! For the real-data-cases only.
1772 REAL , INTENT(OUT) :: p00 , t00 , a , tiso
1773 CALL nl_get_base_pres ( 1 , p00 )
1774 CALL nl_get_base_temp ( 1 , t00 )
1775 CALL nl_get_base_lapse ( 1 , a )
1776 CALL nl_get_iso_temp ( 1 , tiso )
1777 END SUBROUTINE const_module_initialize
1779 !-------------------------------------------------------------------
1781 SUBROUTINE rebalance_driver ( grid )
1783 IMPLICIT NONE
1785 TYPE (domain) :: grid
1787 CALL rebalance( grid &
1788 !
1789 #include "actual_new_args.inc"
1790 !
1791 )
1793 END SUBROUTINE rebalance_driver
1795 !---------------------------------------------------------------------
1797 SUBROUTINE rebalance ( grid &
1798 !
1799 #include "dummy_new_args.inc"
1800 !
1801 )
1802 IMPLICIT NONE
1804 TYPE (domain) :: grid
1806 #include "dummy_new_decl.inc"
1808 TYPE (grid_config_rec_type) :: config_flags
1810 REAL :: p_surf , pd_surf, p_surf_int , pb_int , ht_hold
1811 REAL :: qvf , qvf1 , qvf2
1812 REAL :: p00 , t00 , a , tiso
1813 REAL , DIMENSION(:,:,:) , ALLOCATABLE :: t_init_int
1815 ! Local domain indices and counters.
1817 INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
1819 INTEGER :: &
1820 ids, ide, jds, jde, kds, kde, &
1821 ims, ime, jms, jme, kms, kme, &
1822 its, ite, jts, jte, kts, kte, &
1823 ips, ipe, jps, jpe, kps, kpe, &
1824 i, j, k
1826 REAL :: temp, temp_int
1828 SELECT CASE ( model_data_order )
1829 CASE ( DATA_ORDER_ZXY )
1830 kds = grid%sd31 ; kde = grid%ed31 ;
1831 ids = grid%sd32 ; ide = grid%ed32 ;
1832 jds = grid%sd33 ; jde = grid%ed33 ;
1834 kms = grid%sm31 ; kme = grid%em31 ;
1835 ims = grid%sm32 ; ime = grid%em32 ;
1836 jms = grid%sm33 ; jme = grid%em33 ;
1838 kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch
1839 its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch
1840 jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
1842 CASE ( DATA_ORDER_XYZ )
1843 ids = grid%sd31 ; ide = grid%ed31 ;
1844 jds = grid%sd32 ; jde = grid%ed32 ;
1845 kds = grid%sd33 ; kde = grid%ed33 ;
1847 ims = grid%sm31 ; ime = grid%em31 ;
1848 jms = grid%sm32 ; jme = grid%em32 ;
1849 kms = grid%sm33 ; kme = grid%em33 ;
1851 its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
1852 jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch
1853 kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch
1855 CASE ( DATA_ORDER_XZY )
1856 ids = grid%sd31 ; ide = grid%ed31 ;
1857 kds = grid%sd32 ; kde = grid%ed32 ;
1858 jds = grid%sd33 ; jde = grid%ed33 ;
1860 ims = grid%sm31 ; ime = grid%em31 ;
1861 kms = grid%sm32 ; kme = grid%em32 ;
1862 jms = grid%sm33 ; jme = grid%em33 ;
1864 its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
1865 kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch
1866 jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
1868 END SELECT
1870 ALLOCATE ( t_init_int(ims:ime,kms:kme,jms:jme) )
1872 ! Some of the many weird geopotential initializations that we'll see today: grid%ph0 is total,
1873 ! and grid%ph_2 is a perturbation from the base state geopotential. We set the base geopotential
1874 ! at the lowest level to terrain elevation * gravity.
1876 DO j=jts,jte
1877 DO i=its,ite
1878 grid%ph0(i,1,j) = grid%ht_fine(i,j) * g
1879 grid%ph_2(i,1,j) = 0.
1880 END DO
1881 END DO
1883 ! To define the base state, we call a USER MODIFIED routine to set the three
1884 ! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K),
1885 ! and A (temperature difference, from 1000 mb to 300 mb, K).
1887 CALL const_module_initialize ( p00 , t00 , a , tiso )
1889 ! Base state potential temperature and inverse density (alpha = 1/rho) from
1890 ! the half eta levels and the base-profile surface pressure. Compute 1/rho
1891 ! from equation of state. The potential temperature is a perturbation from t0.
1893 DO j = jts, MIN(jte,jde-1)
1894 DO i = its, MIN(ite,ide-1)
1896 ! Base state pressure is a function of eta level and terrain, only, plus
1897 ! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
1898 ! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
1899 ! The fine grid terrain is ht_fine, the interpolated is grid%ht.
1901 p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht_fine(i,j)/a/r_d ) **0.5 )
1902 p_surf_int = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j) /a/r_d ) **0.5 )
1904 DO k = 1, kte-1
1905 grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
1906 pb_int = grid%znu(k)*(p_surf_int - grid%p_top) + grid%p_top
1907 temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
1908 ! temp = t00 + A*LOG(pb/p00)
1909 grid%t_init(i,k,j) = temp*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
1910 ! 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
1911 temp_int = MAX ( tiso, t00 + A*LOG(pb_int /p00) )
1912 t_init_int(i,k,j)= temp_int*(p00/pb_int )**(r_d/cp) - t0
1913 ! t_init_int(i,k,j)= (t00 + A*LOG(pb_int /p00))*(p00/pb_int )**(r_d/cp) - t0
1914 grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
1915 END DO
1917 ! Base state mu is defined as base state surface pressure minus grid%p_top
1919 grid%mub(i,j) = p_surf - grid%p_top
1921 ! Dry surface pressure is defined as the following (this mu is from the input file
1922 ! computed from the dry pressure). Here the dry pressure is just reconstituted.
1924 pd_surf = ( grid%mub(i,j) + grid%mu_2(i,j) ) + grid%p_top
1926 ! Integrate base geopotential, starting at terrain elevation. This assures that
1927 ! the base state is in exact hydrostatic balance with respect to the model equations.
1928 ! This field is on full levels.
1930 grid%phb(i,1,j) = grid%ht_fine(i,j) * g
1931 DO k = 2,kte
1932 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)
1933 END DO
1934 END DO
1935 END DO
1937 ! Replace interpolated terrain with fine grid values.
1939 DO j = jts, MIN(jte,jde-1)
1940 DO i = its, MIN(ite,ide-1)
1941 grid%ht(i,j) = grid%ht_fine(i,j)
1942 END DO
1943 END DO
1945 ! Perturbation fields.
1947 DO j = jts, min(jde-1,jte)
1948 DO i = its, min(ide-1,ite)
1950 ! The potential temperature is THETAnest = THETAinterp + ( TBARnest - TBARinterp)
1952 DO k = 1 , kde-1
1953 grid%t_2(i,k,j) = grid%t_2(i,k,j) + ( grid%t_init(i,k,j) - t_init_int(i,k,j) )
1954 END DO
1956 ! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
1957 ! equation) down from the top to get the pressure perturbation. First get the pressure
1958 ! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
1960 k = kte-1
1962 qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
1963 qvf2 = 1./(1.+qvf1)
1964 qvf1 = qvf1*qvf2
1966 grid%p(i,k,j) = - 0.5*(grid%mu_2(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
1967 qvf = 1. + rvovrd*moist(i,k,j,P_QV)
1968 grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf* &
1969 (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
1970 grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
1972 ! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
1973 ! inverse density fields (total and perturbation).
1975 DO k=kte-2,1,-1
1976 qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
1977 qvf2 = 1./(1.+qvf1)
1978 qvf1 = qvf1*qvf2
1979 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)
1980 qvf = 1. + rvovrd*moist(i,k,j,P_QV)
1981 grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf* &
1982 (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
1983 grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
1984 END DO
1986 ! This is the hydrostatic equation used in the model after the small timesteps. In
1987 ! the model, grid%al (inverse density) is computed from the geopotential.
1989 DO k = 2,kte
1990 grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
1991 grid%dnw(k-1) * ( (grid%mub(i,j)+grid%mu_2(i,j))*grid%al(i,k-1,j) &
1992 + grid%mu_2(i,j)*grid%alb(i,k-1,j) )
1993 grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
1994 END DO
1996 END DO
1997 END DO
1999 DEALLOCATE ( t_init_int )
2001 ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
2002 #ifdef DM_PARALLEL
2003 # include "HALO_EM_INIT_1.inc"
2004 # include "HALO_EM_INIT_2.inc"
2005 # include "HALO_EM_INIT_3.inc"
2006 # include "HALO_EM_INIT_4.inc"
2007 # include "HALO_EM_INIT_5.inc"
2008 #endif
2009 END SUBROUTINE rebalance
2011 !---------------------------------------------------------------------
2013 RECURSIVE SUBROUTINE find_my_parent ( grid_ptr_in , grid_ptr_out , id_i_am , id_wanted , found_the_id )
2015 USE module_domain
2017 TYPE(domain) , POINTER :: grid_ptr_in , grid_ptr_out
2018 TYPE(domain) , POINTER :: grid_ptr_sibling
2019 INTEGER :: id_wanted , id_i_am
2020 LOGICAL :: found_the_id
2022 found_the_id = .FALSE.
2023 grid_ptr_sibling => grid_ptr_in
2024 DO WHILE ( ASSOCIATED ( grid_ptr_sibling ) )
2026 IF ( grid_ptr_sibling%grid_id .EQ. id_wanted ) THEN
2027 found_the_id = .TRUE.
2028 grid_ptr_out => grid_ptr_sibling
2029 RETURN
2030 ELSE IF ( grid_ptr_sibling%num_nests .GT. 0 ) THEN
2031 grid_ptr_sibling => grid_ptr_sibling%nests(1)%ptr
2032 CALL find_my_parent ( grid_ptr_sibling , grid_ptr_out , id_i_am , id_wanted , found_the_id )
2033 ELSE
2034 grid_ptr_sibling => grid_ptr_sibling%sibling
2035 END IF
2037 END DO
2039 END SUBROUTINE find_my_parent
2041 #endif
2043 !---------------------------------------------------------------------
2045 #ifdef VERT_UNIT
2047 !This is a main program for a small unit test for the vertical interpolation.
2049 program vint
2051 implicit none
2053 integer , parameter :: ij = 3
2054 integer , parameter :: keta = 30
2055 integer , parameter :: kgen =20
2057 integer :: ids , ide , jds , jde , kds , kde , &
2058 ims , ime , jms , jme , kms , kme , &
2059 its , ite , jts , jte , kts , kte
2061 integer :: generic
2063 real , dimension(1:ij,kgen,1:ij) :: fo , po
2064 real , dimension(1:ij,1:keta,1:ij) :: fn_calc , fn_interp , pn
2066 integer, parameter :: interp_type = 1 ! 2
2067 ! integer, parameter :: lagrange_order = 2 ! 1
2068 integer :: lagrange_order
2069 logical, parameter :: lowest_lev_from_sfc = .FALSE. ! .TRUE.
2070 logical, parameter :: use_levels_below_ground = .FALSE. ! .TRUE.
2071 logical, parameter :: use_surface = .FALSE. ! .TRUE.
2072 real , parameter :: zap_close_levels = 500. ! 100.
2073 integer, parameter :: force_sfc_in_vinterp = 0 ! 6
2075 integer :: k
2077 ids = 1 ; ide = ij ; jds = 1 ; jde = ij ; kds = 1 ; kde = keta
2078 ims = 1 ; ime = ij ; jms = 1 ; jme = ij ; kms = 1 ; kme = keta
2079 its = 1 ; ite = ij ; jts = 1 ; jte = ij ; kts = 1 ; kte = keta
2081 generic = kgen
2083 print *,' '
2084 print *,'------------------------------------'
2085 print *,'UNIT TEST FOR VERTICAL INTERPOLATION'
2086 print *,'------------------------------------'
2087 print *,' '
2088 do lagrange_order = 1 , 2
2089 print *,' '
2090 print *,'------------------------------------'
2091 print *,'Lagrange Order = ',lagrange_order
2092 print *,'------------------------------------'
2093 print *,' '
2094 call fillitup ( fo , po , fn_calc , pn , &
2095 ids , ide , jds , jde , kds , kde , &
2096 ims , ime , jms , jme , kms , kme , &
2097 its , ite , jts , jte , kts , kte , &
2098 generic , lagrange_order )
2100 print *,' '
2101 print *,'Level Pressure Field'
2102 print *,' (Pa) (generic)'
2103 print *,'------------------------------------'
2104 print *,' '
2105 do k = 1 , generic
2106 write (*,fmt='(i2,2x,f12.3,1x,g15.8)' ) &
2107 k,po(2,k,2),fo(2,k,2)
2108 end do
2109 print *,' '
2111 call vert_interp ( fo , po , fn_interp , pn , &
2112 generic , 'T' , &
2113 interp_type , lagrange_order , &
2114 lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
2115 zap_close_levels , force_sfc_in_vinterp , &
2116 ids , ide , jds , jde , kds , kde , &
2117 ims , ime , jms , jme , kms , kme , &
2118 its , ite , jts , jte , kts , kte )
2120 print *,'Multi-Order Interpolator'
2121 print *,'------------------------------------'
2122 print *,' '
2123 print *,'Level Pressure Field Field Field'
2124 print *,' (Pa) Calc Interp Diff'
2125 print *,'------------------------------------'
2126 print *,' '
2127 do k = kts , kte-1
2128 write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) &
2129 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)
2130 end do
2132 call vert_interp_old ( fo , po , fn_interp , pn , &
2133 generic , 'T' , &
2134 interp_type , lagrange_order , &
2135 lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
2136 zap_close_levels , force_sfc_in_vinterp , &
2137 ids , ide , jds , jde , kds , kde , &
2138 ims , ime , jms , jme , kms , kme , &
2139 its , ite , jts , jte , kts , kte )
2141 print *,'Linear Interpolator'
2142 print *,'------------------------------------'
2143 print *,' '
2144 print *,'Level Pressure Field Field Field'
2145 print *,' (Pa) Calc Interp Diff'
2146 print *,'------------------------------------'
2147 print *,' '
2148 do k = kts , kte-1
2149 write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) &
2150 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)
2151 end do
2152 end do
2154 end program vint
2156 subroutine wrf_error_fatal (string)
2157 character (len=*) :: string
2158 print *,string
2159 stop
2160 end subroutine wrf_error_fatal
2162 subroutine fillitup ( fo , po , fn , pn , &
2163 ids , ide , jds , jde , kds , kde , &
2164 ims , ime , jms , jme , kms , kme , &
2165 its , ite , jts , jte , kts , kte , &
2166 generic , lagrange_order )
2168 implicit none
2170 integer , intent(in) :: ids , ide , jds , jde , kds , kde , &
2171 ims , ime , jms , jme , kms , kme , &
2172 its , ite , jts , jte , kts , kte
2174 integer , intent(in) :: generic , lagrange_order
2176 real , dimension(ims:ime,generic,jms:jme) , intent(out) :: fo , po
2177 real , dimension(ims:ime,kms:kme,jms:jme) , intent(out) :: fn , pn
2179 integer :: i , j , k
2181 real , parameter :: piov2 = 3.14159265358 / 2.
2183 k = 1
2184 do j = jts , jte
2185 do i = its , ite
2186 po(i,k,j) = 102000.
2187 end do
2188 end do
2190 do k = 2 , generic
2191 do j = jts , jte
2192 do i = its , ite
2193 po(i,k,j) = ( 5000. * ( 1 - (k-1) ) + 100000. * ( (k-1) - (generic-1) ) ) / (1. - real(generic-1) )
2194 end do
2195 end do
2196 end do
2198 if ( lagrange_order .eq. 1 ) then
2199 do k = 1 , generic
2200 do j = jts , jte
2201 do i = its , ite
2202 fo(i,k,j) = po(i,k,j)
2203 ! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. )
2204 end do
2205 end do
2206 end do
2207 else if ( lagrange_order .eq. 2 ) then
2208 do k = 1 , generic
2209 do j = jts , jte
2210 do i = its , ite
2211 fo(i,k,j) = (((po(i,k,j)-5000.)/102000.)*((102000.-po(i,k,j))/102000.))*102000.
2212 ! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. )
2213 end do
2214 end do
2215 end do
2216 end if
2218 !!!!!!!!!!!!
2220 do k = kts , kte
2221 do j = jts , jte
2222 do i = its , ite
2223 pn(i,k,j) = ( 5000. * ( 0 - (k-1) ) + 102000. * ( (k-1) - (kte-1) ) ) / (-1. * real(kte-1) )
2224 end do
2225 end do
2226 end do
2228 do k = kts , kte-1
2229 do j = jts , jte
2230 do i = its , ite
2231 pn(i,k,j) = ( pn(i,k,j) + pn(i,k+1,j) ) /2.
2232 end do
2233 end do
2234 end do
2237 if ( lagrange_order .eq. 1 ) then
2238 do k = kts , kte-1
2239 do j = jts , jte
2240 do i = its , ite
2241 fn(i,k,j) = pn(i,k,j)
2242 ! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. )
2243 end do
2244 end do
2245 end do
2246 else if ( lagrange_order .eq. 2 ) then
2247 do k = kts , kte-1
2248 do j = jts , jte
2249 do i = its , ite
2250 fn(i,k,j) = (((pn(i,k,j)-5000.)/102000.)*((102000.-pn(i,k,j))/102000.))*102000.
2251 ! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. )
2252 end do
2253 end do
2254 end do
2255 end if
2257 end subroutine fillitup
2259 #endif
2261 !---------------------------------------------------------------------
2263 SUBROUTINE vert_interp ( fo , po , fnew , pnu , &
2264 generic , var_type , &
2265 interp_type , lagrange_order , extrap_type , &
2266 lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
2267 zap_close_levels , force_sfc_in_vinterp , &
2268 ids , ide , jds , jde , kds , kde , &
2269 ims , ime , jms , jme , kms , kme , &
2270 its , ite , jts , jte , kts , kte )
2272 ! Vertically interpolate the new field. The original field on the original
2273 ! pressure levels is provided, and the new pressure surfaces to interpolate to.
2275 IMPLICIT NONE
2277 INTEGER , INTENT(IN) :: interp_type , lagrange_order , extrap_type
2278 LOGICAL , INTENT(IN) :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
2279 REAL , INTENT(IN) :: zap_close_levels
2280 INTEGER , INTENT(IN) :: force_sfc_in_vinterp
2281 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2282 ims , ime , jms , jme , kms , kme , &
2283 its , ite , jts , jte , kts , kte
2284 INTEGER , INTENT(IN) :: generic
2286 CHARACTER (LEN=1) :: var_type
2288 REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: fo , po
2289 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu
2290 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew
2292 REAL , DIMENSION(ims:ime,generic,jms:jme) :: forig , porig
2293 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew
2295 ! Local vars
2297 INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2 , knext
2298 INTEGER :: istart , iend , jstart , jend , kstart , kend
2299 INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below
2300 INTEGER , DIMENSION(ims:ime ) :: ks
2301 INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc
2302 INTEGER :: count , zap , zap_below , zap_above , kst , kcount
2303 INTEGER :: kinterp_start , kinterp_end , sfc_level
2305 LOGICAL :: any_below_ground
2307 REAL :: p1 , p2 , pn, hold
2308 REAL , DIMENSION(1:generic) :: ordered_porig , ordered_forig
2309 REAL , DIMENSION(kts:kte) :: ordered_pnew , ordered_fnew
2311 ! Horiontal loop bounds for different variable types.
2313 IF ( var_type .EQ. 'U' ) THEN
2314 istart = its
2315 iend = ite
2316 jstart = jts
2317 jend = MIN(jde-1,jte)
2318 kstart = kts
2319 kend = kte-1
2320 DO j = jstart,jend
2321 DO k = 1,generic
2322 DO i = MAX(ids+1,its) , MIN(ide-1,ite)
2323 porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5
2324 END DO
2325 END DO
2326 IF ( ids .EQ. its ) THEN
2327 DO k = 1,generic
2328 porig(its,k,j) = po(its,k,j)
2329 END DO
2330 END IF
2331 IF ( ide .EQ. ite ) THEN
2332 DO k = 1,generic
2333 porig(ite,k,j) = po(ite-1,k,j)
2334 END DO
2335 END IF
2337 DO k = kstart,kend
2338 DO i = MAX(ids+1,its) , MIN(ide-1,ite)
2339 pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5
2340 END DO
2341 END DO
2342 IF ( ids .EQ. its ) THEN
2343 DO k = kstart,kend
2344 pnew(its,k,j) = pnu(its,k,j)
2345 END DO
2346 END IF
2347 IF ( ide .EQ. ite ) THEN
2348 DO k = kstart,kend
2349 pnew(ite,k,j) = pnu(ite-1,k,j)
2350 END DO
2351 END IF
2352 END DO
2353 ELSE IF ( var_type .EQ. 'V' ) THEN
2354 istart = its
2355 iend = MIN(ide-1,ite)
2356 jstart = jts
2357 jend = jte
2358 kstart = kts
2359 kend = kte-1
2360 DO i = istart,iend
2361 DO k = 1,generic
2362 DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
2363 porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5
2364 END DO
2365 END DO
2366 IF ( jds .EQ. jts ) THEN
2367 DO k = 1,generic
2368 porig(i,k,jts) = po(i,k,jts)
2369 END DO
2370 END IF
2371 IF ( jde .EQ. jte ) THEN
2372 DO k = 1,generic
2373 porig(i,k,jte) = po(i,k,jte-1)
2374 END DO
2375 END IF
2377 DO k = kstart,kend
2378 DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
2379 pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5
2380 END DO
2381 END DO
2382 IF ( jds .EQ. jts ) THEN
2383 DO k = kstart,kend
2384 pnew(i,k,jts) = pnu(i,k,jts)
2385 END DO
2386 END IF
2387 IF ( jde .EQ. jte ) THEN
2388 DO k = kstart,kend
2389 pnew(i,k,jte) = pnu(i,k,jte-1)
2390 END DO
2391 END IF
2392 END DO
2393 ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN
2394 istart = its
2395 iend = MIN(ide-1,ite)
2396 jstart = jts
2397 jend = MIN(jde-1,jte)
2398 kstart = kts
2399 kend = kte
2400 DO j = jstart,jend
2401 DO k = 1,generic
2402 DO i = istart,iend
2403 porig(i,k,j) = po(i,k,j)
2404 END DO
2405 END DO
2407 DO k = kstart,kend
2408 DO i = istart,iend
2409 pnew(i,k,j) = pnu(i,k,j)
2410 END DO
2411 END DO
2412 END DO
2413 ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN
2414 istart = its
2415 iend = MIN(ide-1,ite)
2416 jstart = jts
2417 jend = MIN(jde-1,jte)
2418 kstart = kts
2419 kend = kte-1
2420 DO j = jstart,jend
2421 DO k = 1,generic
2422 DO i = istart,iend
2423 porig(i,k,j) = po(i,k,j)
2424 END DO
2425 END DO
2427 DO k = kstart,kend
2428 DO i = istart,iend
2429 pnew(i,k,j) = pnu(i,k,j)
2430 END DO
2431 END DO
2432 END DO
2433 ELSE
2434 istart = its
2435 iend = MIN(ide-1,ite)
2436 jstart = jts
2437 jend = MIN(jde-1,jte)
2438 kstart = kts
2439 kend = kte-1
2440 DO j = jstart,jend
2441 DO k = 1,generic
2442 DO i = istart,iend
2443 porig(i,k,j) = po(i,k,j)
2444 END DO
2445 END DO
2447 DO k = kstart,kend
2448 DO i = istart,iend
2449 pnew(i,k,j) = pnu(i,k,j)
2450 END DO
2451 END DO
2452 END DO
2453 END IF
2455 DO j = jstart , jend
2457 ! The lowest level is the surface. Levels 2 through "generic" are supposed to
2458 ! be "bottom-up". Flip if they are not. This is based on the input pressure
2459 ! array.
2461 IF ( porig(its,2,j) .LT. porig(its,generic,j) ) THEN
2462 DO kn = 2 , ( generic + 1 ) / 2
2463 DO i = istart , iend
2464 hold = porig(i,kn,j)
2465 porig(i,kn,j) = porig(i,generic+2-kn,j)
2466 porig(i,generic+2-kn,j) = hold
2467 forig(i,kn,j) = fo (i,generic+2-kn,j)
2468 forig(i,generic+2-kn,j) = fo (i,kn,j)
2469 END DO
2470 DO i = istart , iend
2471 forig(i,1,j) = fo (i,1,j)
2472 END DO
2473 END DO
2474 ELSE
2475 DO kn = 1 , generic
2476 DO i = istart , iend
2477 forig(i,kn,j) = fo (i,kn,j)
2478 END DO
2479 END DO
2480 END IF
2482 ! Skip all of the levels below ground in the original data based upon the surface pressure.
2483 ! The ko_above_sfc is the index in the pressure array that is above the surface. If there
2484 ! are no levels underground, this is index = 2. The remaining levels are eligible for use
2485 ! in the vertical interpolation.
2487 DO i = istart , iend
2488 ko_above_sfc(i) = -1
2489 END DO
2490 DO ko = kstart+1 , kend
2491 DO i = istart , iend
2492 IF ( ko_above_sfc(i) .EQ. -1 ) THEN
2493 IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN
2494 ko_above_sfc(i) = ko
2495 END IF
2496 END IF
2497 END DO
2498 END DO
2500 ! Piece together columns of the original input data. Pass the vertical columns to
2501 ! the iterpolator.
2503 DO i = istart , iend
2505 ! If the surface value is in the middle of the array, three steps: 1) do the
2506 ! values below the ground (this is just to catch the occasional value that is
2507 ! inconsistently below the surface based on input data), 2) do the surface level, then
2508 ! 3) add in the levels that are above the surface. For the levels next to the surface,
2509 ! we check to remove any levels that are "too close". When building the column of input
2510 ! pressures, we also attend to the request for forcing the surface analysis to be used
2511 ! in a few lower eta-levels.
2513 ! Fill in the column from up to the level just below the surface with the input
2514 ! presssure and the input field (orig or old, which ever). For an isobaric input
2515 ! file, this data is isobaric.
2517 ! How many levels have we skipped in the input column.
2519 zap = 0
2520 zap_below = 0
2521 zap_above = 0
2523 IF ( ko_above_sfc(i) .GT. 2 ) THEN
2524 count = 1
2525 DO ko = 2 , ko_above_sfc(i)-1
2526 ordered_porig(count) = porig(i,ko,j)
2527 ordered_forig(count) = forig(i,ko,j)
2528 count = count + 1
2529 END DO
2531 ! Make sure the pressure just below the surface is not "too close", this
2532 ! will cause havoc with the higher order interpolators. In case of a "too close"
2533 ! instance, we toss out the offending level (NOT the surface one) by simply
2534 ! decrementing the accumulating loop counter.
2536 IF ( ordered_porig(count-1) - porig(i,1,j) .LT. zap_close_levels ) THEN
2537 count = count -1
2538 zap = 1
2539 zap_below = 1
2540 END IF
2542 ! Add in the surface values.
2544 ordered_porig(count) = porig(i,1,j)
2545 ordered_forig(count) = forig(i,1,j)
2546 count = count + 1
2548 ! A usual way to do the vertical interpolation is to pay more attention to the
2549 ! surface data. Why? Well it has about 20x the density as the upper air, so we
2550 ! hope the analysis is better there. We more strongly use this data by artificially
2551 ! tossing out levels above the surface that are beneath a certain number of prescribed
2552 ! eta levels at this (i,j). The "zap" value is how many levels of input we are
2553 ! removing, which is used to tell the interpolator how many valid values are in
2554 ! the column. The "count" value is the increment to the index of levels, and is
2555 ! only used for assignments.
2557 IF ( force_sfc_in_vinterp .GT. 0 ) THEN
2559 ! Get the pressure at the eta level. We want to remove all input pressure levels
2560 ! between the level above the surface to the pressure at this eta surface. That
2561 ! forces the surface value to be used through the selected eta level. Keep track
2562 ! of two things: the level to use above the eta levels, and how many levels we are
2563 ! skipping.
2565 knext = ko_above_sfc(i)
2566 find_level : DO ko = ko_above_sfc(i) , generic
2567 IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN
2568 knext = ko
2569 exit find_level
2570 ELSE
2571 zap = zap + 1
2572 zap_above = zap_above + 1
2573 END IF
2574 END DO find_level
2576 ! No request for special interpolation, so we just assign the next level to use
2577 ! above the surface as, ta da, the first level above the surface. I know, wow.
2579 ELSE
2580 knext = ko_above_sfc(i)
2581 END IF
2583 ! One more time, make sure the pressure just above the surface is not "too close", this
2584 ! will cause havoc with the higher order interpolators. In case of a "too close"
2585 ! instance, we toss out the offending level above the surface (NOT the surface one) by simply
2586 ! incrementing the loop counter. Here, count-1 is the surface level and knext is either
2587 ! the next level up OR it is the level above the prescribed number of eta surfaces.
2589 IF ( ordered_porig(count-1) - porig(i,knext,j) .LT. zap_close_levels ) THEN
2590 kst = knext+1
2591 zap = zap + 1
2592 zap_above = zap_above + 1
2593 ELSE
2594 kst = knext
2595 END IF
2597 DO ko = kst , generic
2598 ordered_porig(count) = porig(i,ko,j)
2599 ordered_forig(count) = forig(i,ko,j)
2600 count = count + 1
2601 END DO
2603 ! This is easy, the surface is the lowest level, just stick them in, in this order. OK,
2604 ! there are a couple of subtleties. We have to check for that special interpolation that
2605 ! skips some input levels so that the surface is used for the lowest few eta levels. Also,
2606 ! we must macke sure that we still do not have levels that are "too close" together.
2608 ELSE
2610 ! Initialize no input levels have yet been removed from consideration.
2612 zap = 0
2614 ! The surface is the lowest level, so it gets set right away to location 1.
2616 ordered_porig(1) = porig(i,1,j)
2617 ordered_forig(1) = forig(i,1,j)
2619 ! We start filling in the array at loc 2, as in just above the level we just stored.
2621 count = 2
2623 ! Are we forcing the interpolator to skip valid input levels so that the
2624 ! surface data is used through more levels? Essentially as above.
2626 IF ( force_sfc_in_vinterp .GT. 0 ) THEN
2627 knext = 2
2628 find_level2: DO ko = 2 , generic
2629 IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN
2630 knext = ko
2631 exit find_level2
2632 ELSE
2633 zap = zap + 1
2634 zap_above = zap_above + 1
2635 END IF
2636 END DO find_level2
2637 ELSE
2638 knext = 2
2639 END IF
2641 ! Fill in the data above the surface. The "knext" index is either the one
2642 ! just above the surface OR it is the index associated with the level that
2643 ! is just above the pressure at this (i,j) of the top eta level that is to
2644 ! be directly impacted with the surface level in interpolation.
2646 DO ko = knext , generic
2647 IF ( ordered_porig(count-1) - porig(i,ko,j) .LT. zap_close_levels ) THEN
2648 zap = zap + 1
2649 zap_above = zap_above + 1
2650 CYCLE
2651 END IF
2652 ordered_porig(count) = porig(i,ko,j)
2653 ordered_forig(count) = forig(i,ko,j)
2654 count = count + 1
2655 END DO
2657 END IF
2659 ! Now get the column of the "new" pressure data. So, this one is easy.
2661 DO kn = kstart , kend
2662 ordered_pnew(kn) = pnew(i,kn,j)
2663 END DO
2665 ! How many levels (count) are we shipping to the Lagrange interpolator.
2667 IF ( ( use_levels_below_ground ) .AND. ( use_surface ) ) THEN
2669 ! Use all levels, including the input surface, and including the pressure
2670 ! levels below ground. We know to stop when we have reached the top of
2671 ! the input pressure data.
2673 count = 0
2674 find_how_many_1 : DO ko = 1 , generic
2675 IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN
2676 count = count + 1
2677 EXIT find_how_many_1
2678 ELSE
2679 count = count + 1
2680 END IF
2681 END DO find_how_many_1
2682 kinterp_start = 1
2683 kinterp_end = kinterp_start + count - 1
2685 ELSE IF ( ( use_levels_below_ground ) .AND. ( .NOT. use_surface ) ) THEN
2687 ! Use all levels (excluding the input surface) and including the pressure
2688 ! levels below ground. We know to stop when we have reached the top of
2689 ! the input pressure data.
2691 count = 0
2692 find_sfc_2 : DO ko = 1 , generic
2693 IF ( porig(i,1,j) .EQ. ordered_porig(ko) ) THEN
2694 sfc_level = ko
2695 EXIT find_sfc_2
2696 END IF
2697 END DO find_sfc_2
2699 DO ko = sfc_level , generic-1
2700 ordered_porig(ko) = ordered_porig(ko+1)
2701 ordered_forig(ko) = ordered_forig(ko+1)
2702 END DO
2703 ordered_porig(generic) = 1.E-5
2704 ordered_forig(generic) = 1.E10
2706 count = 0
2707 find_how_many_2 : DO ko = 1 , generic
2708 IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN
2709 count = count + 1
2710 EXIT find_how_many_2
2711 ELSE
2712 count = count + 1
2713 END IF
2714 END DO find_how_many_2
2715 kinterp_start = 1
2716 kinterp_end = kinterp_start + count - 1
2718 ELSE IF ( ( .NOT. use_levels_below_ground ) .AND. ( use_surface ) ) THEN
2720 ! Use all levels above the input surface pressure.
2722 kcount = ko_above_sfc(i)-1-zap_below
2723 count = 0
2724 DO ko = 1 , generic
2725 IF ( porig(i,ko,j) .EQ. ordered_porig(kcount) ) THEN
2726 ! write (6,fmt='(f11.3,f11.3,g11.5)') porig(i,ko,j),ordered_porig(kcount),ordered_forig(kcount)
2727 kcount = kcount + 1
2728 count = count + 1
2729 ELSE
2730 ! write (6,fmt='(f11.3 )') porig(i,ko,j)
2731 END IF
2732 END DO
2733 kinterp_start = ko_above_sfc(i)-1-zap_below
2734 kinterp_end = kinterp_start + count - 1
2736 END IF
2738 ! The polynomials are either in pressure or LOG(pressure).
2740 IF ( interp_type .EQ. 1 ) THEN
2741 CALL lagrange_setup ( var_type , &
2742 ordered_porig(kinterp_start:kinterp_end) , &
2743 ordered_forig(kinterp_start:kinterp_end) , &
2744 count , lagrange_order , extrap_type , &
2745 ordered_pnew(kstart:kend) , ordered_fnew , kend-kstart+1 ,i,j)
2746 ELSE
2747 CALL lagrange_setup ( var_type , &
2748 LOG(ordered_porig(kinterp_start:kinterp_end)) , &
2749 ordered_forig(kinterp_start:kinterp_end) , &
2750 count , lagrange_order , extrap_type , &
2751 LOG(ordered_pnew(kstart:kend)) , ordered_fnew , kend-kstart+1 ,i,j)
2752 END IF
2754 ! Save the computed data.
2756 DO kn = kstart , kend
2757 fnew(i,kn,j) = ordered_fnew(kn)
2758 END DO
2760 ! There may have been a request to have the surface data from the input field
2761 ! to be assigned as to the lowest eta level. This assumes thin layers (usually
2762 ! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V).
2764 IF ( lowest_lev_from_sfc ) THEN
2765 fnew(i,1,j) = forig(i,ko_above_sfc(i)-1,j)
2766 END IF
2768 END DO
2770 END DO
2772 END SUBROUTINE vert_interp
2774 !---------------------------------------------------------------------
2776 SUBROUTINE vert_interp_old ( forig , po , fnew , pnu , &
2777 generic , var_type , &
2778 interp_type , lagrange_order , extrap_type , &
2779 lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
2780 zap_close_levels , force_sfc_in_vinterp , &
2781 ids , ide , jds , jde , kds , kde , &
2782 ims , ime , jms , jme , kms , kme , &
2783 its , ite , jts , jte , kts , kte )
2785 ! Vertically interpolate the new field. The original field on the original
2786 ! pressure levels is provided, and the new pressure surfaces to interpolate to.
2788 IMPLICIT NONE
2790 INTEGER , INTENT(IN) :: interp_type , lagrange_order , extrap_type
2791 LOGICAL , INTENT(IN) :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
2792 REAL , INTENT(IN) :: zap_close_levels
2793 INTEGER , INTENT(IN) :: force_sfc_in_vinterp
2794 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2795 ims , ime , jms , jme , kms , kme , &
2796 its , ite , jts , jte , kts , kte
2797 INTEGER , INTENT(IN) :: generic
2799 CHARACTER (LEN=1) :: var_type
2801 REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: forig , po
2802 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu
2803 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew
2805 REAL , DIMENSION(ims:ime,generic,jms:jme) :: porig
2806 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew
2808 ! Local vars
2810 INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2
2811 INTEGER :: istart , iend , jstart , jend , kstart , kend
2812 INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below
2813 INTEGER , DIMENSION(ims:ime ) :: ks
2814 INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc
2816 LOGICAL :: any_below_ground
2818 REAL :: p1 , p2 , pn
2819 integer vert_extrap
2820 vert_extrap = 0
2822 ! Horiontal loop bounds for different variable types.
2824 IF ( var_type .EQ. 'U' ) THEN
2825 istart = its
2826 iend = ite
2827 jstart = jts
2828 jend = MIN(jde-1,jte)
2829 kstart = kts
2830 kend = kte-1
2831 DO j = jstart,jend
2832 DO k = 1,generic
2833 DO i = MAX(ids+1,its) , MIN(ide-1,ite)
2834 porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5
2835 END DO
2836 END DO
2837 IF ( ids .EQ. its ) THEN
2838 DO k = 1,generic
2839 porig(its,k,j) = po(its,k,j)
2840 END DO
2841 END IF
2842 IF ( ide .EQ. ite ) THEN
2843 DO k = 1,generic
2844 porig(ite,k,j) = po(ite-1,k,j)
2845 END DO
2846 END IF
2848 DO k = kstart,kend
2849 DO i = MAX(ids+1,its) , MIN(ide-1,ite)
2850 pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5
2851 END DO
2852 END DO
2853 IF ( ids .EQ. its ) THEN
2854 DO k = kstart,kend
2855 pnew(its,k,j) = pnu(its,k,j)
2856 END DO
2857 END IF
2858 IF ( ide .EQ. ite ) THEN
2859 DO k = kstart,kend
2860 pnew(ite,k,j) = pnu(ite-1,k,j)
2861 END DO
2862 END IF
2863 END DO
2864 ELSE IF ( var_type .EQ. 'V' ) THEN
2865 istart = its
2866 iend = MIN(ide-1,ite)
2867 jstart = jts
2868 jend = jte
2869 kstart = kts
2870 kend = kte-1
2871 DO i = istart,iend
2872 DO k = 1,generic
2873 DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
2874 porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5
2875 END DO
2876 END DO
2877 IF ( jds .EQ. jts ) THEN
2878 DO k = 1,generic
2879 porig(i,k,jts) = po(i,k,jts)
2880 END DO
2881 END IF
2882 IF ( jde .EQ. jte ) THEN
2883 DO k = 1,generic
2884 porig(i,k,jte) = po(i,k,jte-1)
2885 END DO
2886 END IF
2888 DO k = kstart,kend
2889 DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
2890 pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5
2891 END DO
2892 END DO
2893 IF ( jds .EQ. jts ) THEN
2894 DO k = kstart,kend
2895 pnew(i,k,jts) = pnu(i,k,jts)
2896 END DO
2897 END IF
2898 IF ( jde .EQ. jte ) THEN
2899 DO k = kstart,kend
2900 pnew(i,k,jte) = pnu(i,k,jte-1)
2901 END DO
2902 END IF
2903 END DO
2904 ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN
2905 istart = its
2906 iend = MIN(ide-1,ite)
2907 jstart = jts
2908 jend = MIN(jde-1,jte)
2909 kstart = kts
2910 kend = kte
2911 DO j = jstart,jend
2912 DO k = 1,generic
2913 DO i = istart,iend
2914 porig(i,k,j) = po(i,k,j)
2915 END DO
2916 END DO
2918 DO k = kstart,kend
2919 DO i = istart,iend
2920 pnew(i,k,j) = pnu(i,k,j)
2921 END DO
2922 END DO
2923 END DO
2924 ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN
2925 istart = its
2926 iend = MIN(ide-1,ite)
2927 jstart = jts
2928 jend = MIN(jde-1,jte)
2929 kstart = kts
2930 kend = kte-1
2931 DO j = jstart,jend
2932 DO k = 1,generic
2933 DO i = istart,iend
2934 porig(i,k,j) = po(i,k,j)
2935 END DO
2936 END DO
2938 DO k = kstart,kend
2939 DO i = istart,iend
2940 pnew(i,k,j) = pnu(i,k,j)
2941 END DO
2942 END DO
2943 END DO
2944 ELSE
2945 istart = its
2946 iend = MIN(ide-1,ite)
2947 jstart = jts
2948 jend = MIN(jde-1,jte)
2949 kstart = kts
2950 kend = kte-1
2951 DO j = jstart,jend
2952 DO k = 1,generic
2953 DO i = istart,iend
2954 porig(i,k,j) = po(i,k,j)
2955 END DO
2956 END DO
2958 DO k = kstart,kend
2959 DO i = istart,iend
2960 pnew(i,k,j) = pnu(i,k,j)
2961 END DO
2962 END DO
2963 END DO
2964 END IF
2966 DO j = jstart , jend
2968 ! Skip all of the levels below ground in the original data based upon the surface pressure.
2969 ! The ko_above_sfc is the index in the pressure array that is above the surface. If there
2970 ! are no levels underground, this is index = 2. The remaining levels are eligible for use
2971 ! in the vertical interpolation.
2973 DO i = istart , iend
2974 ko_above_sfc(i) = -1
2975 END DO
2976 DO ko = kstart+1 , kend
2977 DO i = istart , iend
2978 IF ( ko_above_sfc(i) .EQ. -1 ) THEN
2979 IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN
2980 ko_above_sfc(i) = ko
2981 END IF
2982 END IF
2983 END DO
2984 END DO
2986 ! Initialize interpolation location. These are the levels in the original pressure
2987 ! data that are physically below and above the targeted new pressure level.
2989 DO kn = kts , kte
2990 DO i = its , ite
2991 k_above(i,kn) = -1
2992 k_below(i,kn) = -2
2993 END DO
2994 END DO
2996 ! Starting location is no lower than previous found location. This is for O(n logn)
2997 ! and not O(n^2), where n is the number of vertical levels to search.
2999 DO i = its , ite
3000 ks(i) = 1
3001 END DO
3003 ! Find trapping layer for interpolation. The kn index runs through all of the "new"
3004 ! levels of data.
3006 DO kn = kstart , kend
3008 DO i = istart , iend
3010 ! For each "new" level (kn), we search to find the trapping levels in the "orig"
3011 ! data. Most of the time, the "new" levels are the eta surfaces, and the "orig"
3012 ! levels are the input pressure levels.
3014 found_trap_above : DO ko = ks(i) , generic-1
3016 ! Because we can have levels in the interpolation that are not valid,
3017 ! let's toss out any candidate orig pressure values that are below ground
3018 ! based on the surface pressure. If the level =1, then this IS the surface
3019 ! level, so we HAVE to keep that one, but maybe not the ones above. If the
3020 ! level (ks) is NOT=1, then we have to just CYCLE our loop to find a legit
3021 ! below-pressure value. If we are not below ground, then we choose two
3022 ! neighboring levels to test whether they surround the new pressure level.
3024 ! The input trapping levels that we are trying is the surface and the first valid
3025 ! level above the surface.
3027 IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .EQ. 1 ) ) THEN
3028 ko_1 = ko
3029 ko_2 = ko_above_sfc(i)
3031 ! The "below" level is underground, cycle until we get to a valid pressure
3032 ! above ground.
3034 ELSE IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .NE. 1 ) ) THEN
3035 CYCLE found_trap_above
3037 ! The "below" level is above the surface, so we are in the clear to test these
3038 ! two levels out.
3040 ELSE
3041 ko_1 = ko
3042 ko_2 = ko+1
3044 END IF
3046 ! The test of the candidate levels: "below" has to have a larger pressure, and
3047 ! "above" has to have a smaller pressure.
3049 ! OK, we found the correct two surrounding levels. The locations are saved for use in the
3050 ! interpolation.
3052 IF ( ( porig(i,ko_1,j) .GE. pnew(i,kn,j) ) .AND. &
3053 ( porig(i,ko_2,j) .LT. pnew(i,kn,j) ) ) THEN
3054 k_above(i,kn) = ko_2
3055 k_below(i,kn) = ko_1
3056 ks(i) = ko_1
3057 EXIT found_trap_above
3059 ! What do we do is we need to extrapolate the data underground? This happens when the
3060 ! lowest pressure that we have is physically "above" the new target pressure. Our
3061 ! actions depend on the type of variable we are interpolating.
3063 ELSE IF ( porig(i,1,j) .LT. pnew(i,kn,j) ) THEN
3065 ! For horizontal winds and moisture, we keep a constant value under ground.
3067 IF ( ( var_type .EQ. 'U' ) .OR. &
3068 ( var_type .EQ. 'V' ) .OR. &
3069 ( var_type .EQ. 'Q' ) ) THEN
3070 k_above(i,kn) = 1
3071 ks(i) = 1
3073 ! For temperature and height, we extrapolate the data. Hopefully, we are not
3074 ! extrapolating too far. For pressure level input, the eta levels are always
3075 ! contained within the surface to p_top levels, so no extrapolation is ever
3076 ! required.
3078 ELSE IF ( ( var_type .EQ. 'Z' ) .OR. &
3079 ( var_type .EQ. 'T' ) ) THEN
3080 k_above(i,kn) = ko_above_sfc(i)
3081 k_below(i,kn) = 1
3082 ks(i) = 1
3084 ! Just a catch all right now.
3086 ELSE
3087 k_above(i,kn) = 1
3088 ks(i) = 1
3089 END IF
3091 EXIT found_trap_above
3093 ! The other extrapolation that might be required is when we are going above the
3094 ! top level of the input data. Usually this means we chose a P_PTOP value that
3095 ! was inappropriate, and we should stop and let someone fix this mess.
3097 ELSE IF ( porig(i,generic,j) .GT. pnew(i,kn,j) ) THEN
3098 print *,'data is too high, try a lower p_top'
3099 print *,'pnew=',pnew(i,kn,j)
3100 print *,'porig=',porig(i,:,j)
3101 CALL wrf_error_fatal ('requested p_top is higher than input data, lower p_top')
3103 END IF
3104 END DO found_trap_above
3105 END DO
3106 END DO
3108 ! Linear vertical interpolation.
3110 DO kn = kstart , kend
3111 DO i = istart , iend
3112 IF ( k_above(i,kn) .EQ. 1 ) THEN
3113 fnew(i,kn,j) = forig(i,1,j)
3114 ELSE
3115 k2 = MAX ( k_above(i,kn) , 2)
3116 k1 = MAX ( k_below(i,kn) , 1)
3117 IF ( k1 .EQ. k2 ) THEN
3118 CALL wrf_error_fatal ( 'identical values in the interp, bad for divisions' )
3119 END IF
3120 IF ( interp_type .EQ. 1 ) THEN
3121 p1 = porig(i,k1,j)
3122 p2 = porig(i,k2,j)
3123 pn = pnew(i,kn,j)
3124 ELSE IF ( interp_type .EQ. 2 ) THEN
3125 p1 = ALOG(porig(i,k1,j))
3126 p2 = ALOG(porig(i,k2,j))
3127 pn = ALOG(pnew(i,kn,j))
3128 END IF
3129 IF ( ( p1-pn) * (p2-pn) > 0. ) THEN
3130 ! CALL wrf_error_fatal ( 'both trapping pressures are on the same side of the new pressure' )
3131 ! CALL wrf_debug ( 0 , 'both trapping pressures are on the same side of the new pressure' )
3132 vert_extrap = vert_extrap + 1
3133 END IF
3134 fnew(i,kn,j) = ( forig(i,k1,j) * ( p2 - pn ) + &
3135 forig(i,k2,j) * ( pn - p1 ) ) / &
3136 ( p2 - p1 )
3137 END IF
3138 END DO
3139 END DO
3141 search_below_ground : DO kn = kstart , kend
3142 any_below_ground = .FALSE.
3143 DO i = istart , iend
3144 IF ( k_above(i,kn) .EQ. 1 ) THEN
3145 fnew(i,kn,j) = forig(i,1,j)
3146 any_below_ground = .TRUE.
3147 END IF
3148 END DO
3149 IF ( .NOT. any_below_ground ) THEN
3150 EXIT search_below_ground
3151 END IF
3152 END DO search_below_ground
3154 ! There may have been a request to have the surface data from the input field
3155 ! to be assigned as to the lowest eta level. This assumes thin layers (usually
3156 ! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V).
3158 DO i = istart , iend
3159 IF ( lowest_lev_from_sfc ) THEN
3160 fnew(i,1,j) = forig(i,ko_above_sfc(i),j)
3161 END IF
3162 END DO
3164 END DO
3165 print *,'VERT EXTRAP = ', vert_extrap
3167 END SUBROUTINE vert_interp_old
3169 !---------------------------------------------------------------------
3171 SUBROUTINE lagrange_setup ( var_type , all_x , all_y , all_dim , n , extrap_type , &
3172 target_x , target_y , target_dim ,i,j)
3174 ! We call a Lagrange polynomial interpolator. The parallel concerns are put off as this
3175 ! is initially set up for vertical use. The purpose is an input column of pressure (all_x),
3176 ! and the associated pressure level data (all_y). These are assumed to be sorted (ascending
3177 ! or descending, no matter). The locations to be interpolated to are the pressures in
3178 ! target_x, probably the new vertical coordinate values. The field that is output is the
3179 ! target_y, which is defined at the target_x location. Mostly we expect to be 2nd order
3180 ! overlapping polynomials, with only a single 2nd order method near the top and bottom.
3181 ! When n=1, this is linear; when n=2, this is a second order interpolator.
3183 IMPLICIT NONE
3185 CHARACTER (LEN=1) :: var_type
3186 INTEGER , INTENT(IN) :: all_dim , n , extrap_type , target_dim
3187 REAL, DIMENSION(all_dim) , INTENT(IN) :: all_x , all_y
3188 REAL , DIMENSION(target_dim) , INTENT(IN) :: target_x
3189 REAL , DIMENSION(target_dim) , INTENT(OUT) :: target_y
3191 ! Brought in for debug purposes, all of the computations are in a single column.
3193 INTEGER , INTENT(IN) :: i,j
3195 ! Local vars
3197 REAL , DIMENSION(n+1) :: x , y
3198 REAL :: a , b
3199 REAL :: target_y_1 , target_y_2
3200 LOGICAL :: found_loc
3201 INTEGER :: loop , loc_center_left , loc_center_right , ist , iend , target_loop
3202 INTEGER :: vboundb , vboundt
3204 ! Local vars for the problem of extrapolating theta below ground.
3206 REAL :: temp_1 , temp_2 , temp_3 , temp_y
3207 REAL :: depth_of_extrap_in_p , avg_of_extrap_p , temp_extrap_starting_point , dhdp , dh , dt
3208 REAL , PARAMETER :: RovCp = 287. / 1004.
3209 REAL , PARAMETER :: CRC_const1 = 11880.516 ! m
3210 REAL , PARAMETER :: CRC_const2 = 0.1902632 !
3211 REAL , PARAMETER :: CRC_const3 = 0.0065 ! K/km
3213 IF ( all_dim .LT. n+1 ) THEN
3214 print *,'all_dim = ',all_dim
3215 print *,'order = ',n
3216 print *,'i,j = ',i,j
3217 print *,'p array = ',all_x
3218 print *,'f array = ',all_y
3219 print *,'p target= ',target_x
3220 CALL wrf_error_fatal ( 'troubles, the interpolating order is too large for this few input values' )
3221 END IF
3223 IF ( n .LT. 1 ) THEN
3224 CALL wrf_error_fatal ( 'pal, linear is about as low as we go' )
3225 END IF
3227 ! We can pinch in the area of the higher order interpolation with vbound. If
3228 ! vbound = 0, no pinching. If vbound = m, then we make the lower "m" and upper
3229 ! "m" eta levels use a linear interpolation.
3231 vboundb = 4
3232 vboundt = 0
3234 ! Loop over the list of target x and y values.
3236 DO target_loop = 1 , target_dim
3238 ! Find the two trapping x values, and keep the indices.
3240 found_loc = .FALSE.
3241 find_trap : DO loop = 1 , all_dim -1
3242 a = target_x(target_loop) - all_x(loop)
3243 b = target_x(target_loop) - all_x(loop+1)
3244 IF ( a*b .LE. 0.0 ) THEN
3245 loc_center_left = loop
3246 loc_center_right = loop+1
3247 found_loc = .TRUE.
3248 EXIT find_trap
3249 END IF
3250 END DO find_trap
3252 IF ( ( .NOT. found_loc ) .AND. ( target_x(target_loop) .GT. all_x(1) ) ) THEN
3254 ! Isothermal extrapolation.
3256 IF ( ( extrap_type .EQ. 1 ) .AND. ( var_type .EQ. 'T' ) ) THEN
3258 temp_1 = all_y(1) * ( all_x(1) / 100000. ) ** RovCp
3259 target_y(target_loop) = temp_1 * ( 100000. / target_x(target_loop) ) ** RovCp
3261 ! Standard atmosphere -6.5 K/km lapse rate for the extrapolation.
3263 ELSE IF ( ( extrap_type .EQ. 2 ) .AND. ( var_type .EQ. 'T' ) ) THEN
3265 depth_of_extrap_in_p = target_x(target_loop) - all_x(1)
3266 avg_of_extrap_p = ( target_x(target_loop) + all_x(1) ) * 0.5
3267 temp_extrap_starting_point = all_y(1) * ( all_x(1) / 100000. ) ** RovCp
3268 dhdp = CRC_const1 * CRC_const2 * ( avg_of_extrap_p / 100. ) ** ( CRC_const2 - 1. )
3269 dh = dhdp * ( depth_of_extrap_in_p / 100. )
3270 dt = dh * CRC_const3
3271 target_y(target_loop) = ( temp_extrap_starting_point + dt ) * ( 100000. / target_x(target_loop) ) ** RovCp
3273 ! Adiabatic extrapolation for theta.
3275 ELSE IF ( ( extrap_type .EQ. 3 ) .AND. ( var_type .EQ. 'T' ) ) THEN
3277 target_y(target_loop) = all_y(1)
3280 ! Wild extrapolation for non-temperature vars.
3282 ELSE IF ( extrap_type .EQ. 1 ) THEN
3284 target_y(target_loop) = ( all_y(2) * ( target_x(target_loop) - all_x(3) ) + &
3285 all_y(3) * ( all_x(2) - target_x(target_loop) ) ) / &
3286 ( all_x(2) - all_x(3) )
3288 ! Use a constant value below ground.
3290 ELSE IF ( extrap_type .EQ. 2 ) THEN
3292 target_y(target_loop) = all_y(1)
3294 ELSE IF ( extrap_type .EQ. 3 ) THEN
3295 CALL wrf_error_fatal ( 'You are not allowed to use extrap_option #3 for any var except for theta.' )
3297 END IF
3298 CYCLE
3299 ELSE IF ( .NOT. found_loc ) THEN
3300 print *,'i,j = ',i,j
3301 print *,'target pressure and value = ',target_x(target_loop),target_y(target_loop)
3302 DO loop = 1 , all_dim
3303 print *,'column of pressure and value = ',all_x(loop),all_y(loop)
3304 END DO
3305 CALL wrf_error_fatal ( 'troubles, could not find trapping x locations' )
3306 END IF
3308 ! Even or odd order? We can put the value in the middle if this is
3309 ! an odd order interpolator. For the even guys, we'll do it twice
3310 ! and shift the range one index, then get an average.
3312 IF ( MOD(n,2) .NE. 0 ) THEN
3313 IF ( ( loc_center_left -(((n+1)/2)-1) .GE. 1 ) .AND. &
3314 ( loc_center_right+(((n+1)/2)-1) .LE. all_dim ) ) THEN
3315 ist = loc_center_left -(((n+1)/2)-1)
3316 iend = ist + n
3317 CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
3318 ELSE
3319 IF ( .NOT. found_loc ) THEN
3320 CALL wrf_error_fatal ( 'I doubt this will happen, I will only do 2nd order for now' )
3321 END IF
3322 END IF
3324 ELSE IF ( ( MOD(n,2) .EQ. 0 ) .AND. &
3325 ( ( target_loop .GE. 1 + vboundb ) .AND. ( target_loop .LE. target_dim - vboundt ) ) ) THEN
3326 IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. &
3327 ( loc_center_right+(((n )/2) ) .LE. all_dim ) .AND. &
3328 ( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. &
3329 ( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN
3330 ist = loc_center_left -(((n )/2)-1)
3331 iend = ist + n
3332 CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_1 )
3333 ist = loc_center_left -(((n )/2) )
3334 iend = ist + n
3335 CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_2 )
3336 target_y(target_loop) = ( target_y_1 + target_y_2 ) * 0.5
3338 ELSE IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. &
3339 ( loc_center_right+(((n )/2) ) .LE. all_dim ) ) THEN
3340 ist = loc_center_left -(((n )/2)-1)
3341 iend = ist + n
3342 CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
3343 ELSE IF ( ( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. &
3344 ( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN
3345 ist = loc_center_left -(((n )/2) )
3346 iend = ist + n
3347 CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
3348 ELSE
3349 CALL wrf_error_fatal ( 'unauthorized area, you should not be here' )
3350 END IF
3352 ELSE IF ( MOD(n,2) .EQ. 0 ) THEN
3353 ist = loc_center_left
3354 iend = loc_center_right
3355 CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , 1 , target_x(target_loop) , target_y(target_loop) )
3357 END IF
3359 END DO
3361 END SUBROUTINE lagrange_setup
3363 !---------------------------------------------------------------------
3365 SUBROUTINE lagrange_interp ( x , y , n , target_x , target_y )
3367 ! Interpolation using Lagrange polynomials.
3368 ! P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
3369 ! where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
3370 ! ---------------------------------------------
3371 ! (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
3373 IMPLICIT NONE
3375 INTEGER , INTENT(IN) :: n
3376 REAL , DIMENSION(0:n) , INTENT(IN) :: x , y
3377 REAL , INTENT(IN) :: target_x
3379 REAL , INTENT(OUT) :: target_y
3381 ! Local vars
3383 INTEGER :: i , k
3384 REAL :: numer , denom , Px
3385 REAL , DIMENSION(0:n) :: Ln
3387 Px = 0.
3388 DO i = 0 , n
3389 numer = 1.
3390 denom = 1.
3391 DO k = 0 , n
3392 IF ( k .EQ. i ) CYCLE
3393 numer = numer * ( target_x - x(k) )
3394 denom = denom * ( x(i) - x(k) )
3395 END DO
3396 Ln(i) = y(i) * numer / denom
3397 Px = Px + Ln(i)
3398 END DO
3399 target_y = Px
3401 END SUBROUTINE lagrange_interp
3403 #ifndef VERT_UNIT
3404 !---------------------------------------------------------------------
3406 SUBROUTINE p_dry ( mu0 , eta , pdht , pdry , full_levs , &
3407 ids , ide , jds , jde , kds , kde , &
3408 ims , ime , jms , jme , kms , kme , &
3409 its , ite , jts , jte , kts , kte )
3411 ! Compute reference pressure and the reference mu.
3413 IMPLICIT NONE
3415 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3416 ims , ime , jms , jme , kms , kme , &
3417 its , ite , jts , jte , kts , kte
3419 LOGICAL :: full_levs
3421 REAL , DIMENSION(ims:ime, jms:jme) , INTENT(IN) :: mu0
3422 REAL , DIMENSION( kms:kme ) , INTENT(IN) :: eta
3423 REAL :: pdht
3424 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pdry
3426 ! Local vars
3428 INTEGER :: i , j , k
3429 REAL , DIMENSION( kms:kme ) :: eta_h
3431 IF ( full_levs ) THEN
3432 DO j = jts , MIN ( jde-1 , jte )
3433 DO k = kts , kte
3434 DO i = its , MIN (ide-1 , ite )
3435 pdry(i,k,j) = eta(k) * mu0(i,j) + pdht
3436 END DO
3437 END DO
3438 END DO
3440 ELSE
3441 DO k = kts , kte-1
3442 eta_h(k) = ( eta(k) + eta(k+1) ) * 0.5
3443 END DO
3445 DO j = jts , MIN ( jde-1 , jte )
3446 DO k = kts , kte-1
3447 DO i = its , MIN (ide-1 , ite )
3448 pdry(i,k,j) = eta_h(k) * mu0(i,j) + pdht
3449 END DO
3450 END DO
3451 END DO
3452 END IF
3454 END SUBROUTINE p_dry
3456 !---------------------------------------------------------------------
3458 SUBROUTINE p_dts ( pdts , intq , psfc , p_top , &
3459 ids , ide , jds , jde , kds , kde , &
3460 ims , ime , jms , jme , kms , kme , &
3461 its , ite , jts , jte , kts , kte )
3463 ! Compute difference between the dry, total surface pressure and the top pressure.
3465 IMPLICIT NONE
3467 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3468 ims , ime , jms , jme , kms , kme , &
3469 its , ite , jts , jte , kts , kte
3471 REAL , INTENT(IN) :: p_top
3472 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: psfc
3473 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: intq
3474 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(OUT) :: pdts
3476 ! Local vars
3478 INTEGER :: i , j , k
3480 DO j = jts , MIN ( jde-1 , jte )
3481 DO i = its , MIN (ide-1 , ite )
3482 pdts(i,j) = psfc(i,j) - intq(i,j) - p_top
3483 END DO
3484 END DO
3486 END SUBROUTINE p_dts
3488 !---------------------------------------------------------------------
3490 SUBROUTINE p_dhs ( pdhs , ht , p0 , t0 , a , &
3491 ids , ide , jds , jde , kds , kde , &
3492 ims , ime , jms , jme , kms , kme , &
3493 its , ite , jts , jte , kts , kte )
3495 ! Compute dry, hydrostatic surface pressure.
3497 IMPLICIT NONE
3499 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3500 ims , ime , jms , jme , kms , kme , &
3501 its , ite , jts , jte , kts , kte
3503 REAL , DIMENSION(ims:ime, jms:jme) , INTENT(IN) :: ht
3504 REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: pdhs
3506 REAL , INTENT(IN) :: p0 , t0 , a
3508 ! Local vars
3510 INTEGER :: i , j , k
3512 REAL , PARAMETER :: Rd = 287.
3513 REAL , PARAMETER :: g = 9.8
3515 DO j = jts , MIN ( jde-1 , jte )
3516 DO i = its , MIN (ide-1 , ite )
3517 pdhs(i,j) = p0 * EXP ( -t0/a + SQRT ( (t0/a)**2 - 2. * g * ht(i,j)/(a * Rd) ) )
3518 END DO
3519 END DO
3521 END SUBROUTINE p_dhs
3523 !---------------------------------------------------------------------
3525 SUBROUTINE find_p_top ( p , p_top , &
3526 ids , ide , jds , jde , kds , kde , &
3527 ims , ime , jms , jme , kms , kme , &
3528 its , ite , jts , jte , kts , kte )
3530 ! Find the largest pressure in the top level. This is our p_top. We are
3531 ! assuming that the top level is the location where the pressure is a minimum
3532 ! for each column. In cases where the top surface is not isobaric, a
3533 ! communicated value must be shared in the calling routine. Also in cases
3534 ! where the top surface is not isobaric, care must be taken that the new
3535 ! maximum pressure is not greater than the previous value. This test is
3536 ! also handled in the calling routine.
3538 IMPLICIT NONE
3540 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3541 ims , ime , jms , jme , kms , kme , &
3542 its , ite , jts , jte , kts , kte
3544 REAL :: p_top
3545 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
3547 ! Local vars
3549 INTEGER :: i , j , k, min_lev
3551 i = its
3552 j = jts
3553 p_top = p(i,2,j)
3554 min_lev = 2
3555 DO k = 2 , kte
3556 IF ( p_top .GT. p(i,k,j) ) THEN
3557 p_top = p(i,k,j)
3558 min_lev = k
3559 END IF
3560 END DO
3562 k = min_lev
3563 p_top = p(its,k,jts)
3564 DO j = jts , MIN ( jde-1 , jte )
3565 DO i = its , MIN (ide-1 , ite )
3566 p_top = MAX ( p_top , p(i,k,j) )
3567 END DO
3568 END DO
3570 END SUBROUTINE find_p_top
3572 !---------------------------------------------------------------------
3574 SUBROUTINE t_to_theta ( t , p , p00 , &
3575 ids , ide , jds , jde , kds , kde , &
3576 ims , ime , jms , jme , kms , kme , &
3577 its , ite , jts , jte , kts , kte )
3579 ! Compute dry, hydrostatic surface pressure.
3581 IMPLICIT NONE
3583 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3584 ims , ime , jms , jme , kms , kme , &
3585 its , ite , jts , jte , kts , kte
3587 REAL , INTENT(IN) :: p00
3588 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
3589 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: t
3591 ! Local vars
3593 INTEGER :: i , j , k
3595 REAL , PARAMETER :: Rd = 287.
3596 REAL , PARAMETER :: Cp = 1004.
3598 DO j = jts , MIN ( jde-1 , jte )
3599 DO k = kts , kte
3600 DO i = its , MIN (ide-1 , ite )
3601 t(i,k,j) = t(i,k,j) * ( p00 / p(i,k,j) ) ** (Rd / Cp)
3602 END DO
3603 END DO
3604 END DO
3606 END SUBROUTINE t_to_theta
3608 !---------------------------------------------------------------------
3610 SUBROUTINE integ_moist ( q_in , p_in , pd_out , t_in , ght_in , intq , &
3611 ids , ide , jds , jde , kds , kde , &
3612 ims , ime , jms , jme , kms , kme , &
3613 its , ite , jts , jte , kts , kte )
3615 ! Integrate the moisture field vertically. Mostly used to get the total
3616 ! vapor pressure, which can be subtracted from the total pressure to get
3617 ! the dry pressure.
3619 IMPLICIT NONE
3621 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3622 ims , ime , jms , jme , kms , kme , &
3623 its , ite , jts , jte , kts , kte
3625 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: q_in , p_in , t_in , ght_in
3626 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pd_out
3627 REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: intq
3629 ! Local vars
3631 INTEGER :: i , j , k
3632 INTEGER , DIMENSION(ims:ime) :: level_above_sfc
3633 REAL , DIMENSION(ims:ime,jms:jme) :: psfc , tsfc , qsfc, zsfc
3634 REAL , DIMENSION(ims:ime,kms:kme) :: q , p , t , ght, pd
3636 REAL :: rhobar , qbar , dz
3637 REAL :: p1 , p2 , t1 , t2 , q1 , q2 , z1, z2
3639 LOGICAL :: upside_down
3641 REAL , PARAMETER :: Rd = 287.
3642 REAL , PARAMETER :: g = 9.8
3644 ! Get a surface value, always the first level of a 3d field.
3646 DO j = jts , MIN ( jde-1 , jte )
3647 DO i = its , MIN (ide-1 , ite )
3648 psfc(i,j) = p_in(i,kts,j)
3649 tsfc(i,j) = t_in(i,kts,j)
3650 qsfc(i,j) = q_in(i,kts,j)
3651 zsfc(i,j) = ght_in(i,kts,j)
3652 END DO
3653 END DO
3655 IF ( p_in(its,kts+1,jts) .LT. p_in(its,kte,jts) ) THEN
3656 upside_down = .TRUE.
3657 ELSE
3658 upside_down = .FALSE.
3659 END IF
3661 DO j = jts , MIN ( jde-1 , jte )
3663 ! Initialize the integrated quantity of moisture to zero.
3665 DO i = its , MIN (ide-1 , ite )
3666 intq(i,j) = 0.
3667 END DO
3669 IF ( upside_down ) THEN
3670 DO i = its , MIN (ide-1 , ite )
3671 p(i,kts) = p_in(i,kts,j)
3672 t(i,kts) = t_in(i,kts,j)
3673 q(i,kts) = q_in(i,kts,j)
3674 ght(i,kts) = ght_in(i,kts,j)
3675 DO k = kts+1,kte
3676 p(i,k) = p_in(i,kte+2-k,j)
3677 t(i,k) = t_in(i,kte+2-k,j)
3678 q(i,k) = q_in(i,kte+2-k,j)
3679 ght(i,k) = ght_in(i,kte+2-k,j)
3680 END DO
3681 END DO
3682 ELSE
3683 DO i = its , MIN (ide-1 , ite )
3684 DO k = kts,kte
3685 p(i,k) = p_in(i,k ,j)
3686 t(i,k) = t_in(i,k ,j)
3687 q(i,k) = q_in(i,k ,j)
3688 ght(i,k) = ght_in(i,k ,j)
3689 END DO
3690 END DO
3691 END IF
3693 ! Find the first level above the ground. If all of the levels are above ground, such as
3694 ! a terrain following lower coordinate, then the first level above ground is index #2.
3696 DO i = its , MIN (ide-1 , ite )
3697 level_above_sfc(i) = -1
3698 IF ( p(i,kts+1) .LT. psfc(i,j) ) THEN
3699 level_above_sfc(i) = kts+1
3700 ELSE
3701 find_k : DO k = kts+1,kte-1
3702 IF ( ( p(i,k )-psfc(i,j) .GE. 0. ) .AND. &
3703 ( p(i,k+1)-psfc(i,j) .LT. 0. ) ) THEN
3704 level_above_sfc(i) = k+1
3705 EXIT find_k
3706 END IF
3707 END DO find_k
3708 IF ( level_above_sfc(i) .EQ. -1 ) THEN
3709 print *,'i,j = ',i,j
3710 print *,'p = ',p(i,:)
3711 print *,'p sfc = ',psfc(i,j)
3712 CALL wrf_error_fatal ( 'Could not find level above ground')
3713 END IF
3714 END IF
3715 END DO
3717 DO i = its , MIN (ide-1 , ite )
3719 ! Account for the moisture above the ground.
3721 pd(i,kte) = p(i,kte)
3722 DO k = kte-1,level_above_sfc(i),-1
3723 rhobar = ( p(i,k ) / ( Rd * t(i,k ) ) + &
3724 p(i,k+1) / ( Rd * t(i,k+1) ) ) * 0.5
3725 qbar = ( q(i,k ) + q(i,k+1) ) * 0.5
3726 dz = ght(i,k+1) - ght(i,k)
3727 intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz
3728 pd(i,k) = p(i,k) - intq(i,j)
3729 END DO
3731 ! Account for the moisture between the surface and the first level up.
3733 IF ( ( p(i,level_above_sfc(i)-1)-psfc(i,j) .GE. 0. ) .AND. &
3734 ( p(i,level_above_sfc(i) )-psfc(i,j) .LT. 0. ) .AND. &
3735 ( level_above_sfc(i) .GT. kts ) ) THEN
3736 p1 = psfc(i,j)
3737 p2 = p(i,level_above_sfc(i))
3738 t1 = tsfc(i,j)
3739 t2 = t(i,level_above_sfc(i))
3740 q1 = qsfc(i,j)
3741 q2 = q(i,level_above_sfc(i))
3742 z1 = zsfc(i,j)
3743 z2 = ght(i,level_above_sfc(i))
3744 rhobar = ( p1 / ( Rd * t1 ) + &
3745 p2 / ( Rd * t2 ) ) * 0.5
3746 qbar = ( q1 + q2 ) * 0.5
3747 dz = z2 - z1
3748 IF ( dz .GT. 0.1 ) THEN
3749 intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz
3750 END IF
3752 ! Fix the underground values.
3754 DO k = level_above_sfc(i)-1,kts+1,-1
3755 pd(i,k) = p(i,k) - intq(i,j)
3756 END DO
3757 END IF
3758 pd(i,kts) = psfc(i,j) - intq(i,j)
3760 END DO
3762 IF ( upside_down ) THEN
3763 DO i = its , MIN (ide-1 , ite )
3764 pd_out(i,kts,j) = pd(i,kts)
3765 DO k = kts+1,kte
3766 pd_out(i,kte+2-k,j) = pd(i,k)
3767 END DO
3768 END DO
3769 ELSE
3770 DO i = its , MIN (ide-1 , ite )
3771 DO k = kts,kte
3772 pd_out(i,k,j) = pd(i,k)
3773 END DO
3774 END DO
3775 END IF
3777 END DO
3779 END SUBROUTINE integ_moist
3781 !---------------------------------------------------------------------
3783 SUBROUTINE rh_to_mxrat (rh, t, p, q , wrt_liquid , &
3784 ids , ide , jds , jde , kds , kde , &
3785 ims , ime , jms , jme , kms , kme , &
3786 its , ite , jts , jte , kts , kte )
3788 IMPLICIT NONE
3790 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3791 ims , ime , jms , jme , kms , kme , &
3792 its , ite , jts , jte , kts , kte
3794 LOGICAL , INTENT(IN) :: wrt_liquid
3796 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
3797 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
3798 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: q
3800 ! Local vars
3802 INTEGER :: i , j , k
3804 REAL :: ew , q1 , t1
3806 REAL, PARAMETER :: T_REF = 0.0
3807 REAL, PARAMETER :: MW_AIR = 28.966
3808 REAL, PARAMETER :: MW_VAP = 18.0152
3810 REAL, PARAMETER :: A0 = 6.107799961
3811 REAL, PARAMETER :: A1 = 4.436518521e-01
3812 REAL, PARAMETER :: A2 = 1.428945805e-02
3813 REAL, PARAMETER :: A3 = 2.650648471e-04
3814 REAL, PARAMETER :: A4 = 3.031240396e-06
3815 REAL, PARAMETER :: A5 = 2.034080948e-08
3816 REAL, PARAMETER :: A6 = 6.136820929e-11
3818 REAL, PARAMETER :: ES0 = 6.1121
3820 REAL, PARAMETER :: C1 = 9.09718
3821 REAL, PARAMETER :: C2 = 3.56654
3822 REAL, PARAMETER :: C3 = 0.876793
3823 REAL, PARAMETER :: EIS = 6.1071
3824 REAL :: RHS
3825 REAL, PARAMETER :: TF = 273.16
3826 REAL :: TK
3828 REAL :: ES
3829 REAL :: QS
3830 REAL, PARAMETER :: EPS = 0.622
3831 REAL, PARAMETER :: SVP1 = 0.6112
3832 REAL, PARAMETER :: SVP2 = 17.67
3833 REAL, PARAMETER :: SVP3 = 29.65
3834 REAL, PARAMETER :: SVPT0 = 273.15
3836 ! This subroutine computes mixing ratio (q, kg/kg) from basic variables
3837 ! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%).
3838 ! The reference temperature (t_ref, C) is used to describe the temperature
3839 ! at which the liquid and ice phase change occurs.
3841 DO j = jts , MIN ( jde-1 , jte )
3842 DO k = kts , kte
3843 DO i = its , MIN (ide-1 , ite )
3844 rh(i,k,j) = MIN ( MAX ( rh(i,k,j) , 0. ) , 100. )
3845 END DO
3846 END DO
3847 END DO
3849 IF ( wrt_liquid ) THEN
3850 DO j = jts , MIN ( jde-1 , jte )
3851 DO k = kts , kte
3852 DO i = its , MIN (ide-1 , ite )
3854 ! es is reduced by RH here to avoid problems in low-pressure cases
3856 es=.01*rh(i,k,j)*svp1*10.*EXP(svp2*(t(i,k,j)-svpt0)/(t(i,k,j)-svp3))
3857 IF (es .ge. p(i,k,j)/100.)THEN
3858 q(i,k,j)=1.0
3859 print *,'warning: vapor pressure exceeds total pressure '
3860 print *,'setting mixing ratio to 1'
3861 ELSE
3862 q(i,k,j)=eps*es/(p(i,k,j)/100.-es)
3863 ENDIF
3864 END DO
3865 END DO
3866 END DO
3868 ELSE
3869 DO j = jts , MIN ( jde-1 , jte )
3870 DO k = kts , kte
3871 DO i = its , MIN (ide-1 , ite )
3873 t1 = t(i,k,j) - 273.16
3875 ! Obviously dry.
3877 IF ( t1 .lt. -200. ) THEN
3878 q(i,k,j) = 0
3880 ELSE
3882 ! First compute the ambient vapor pressure of water
3884 IF ( ( t1 .GE. t_ref ) .AND. ( t1 .GE. -47.) ) THEN ! liq phase ESLO
3885 ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6)))))
3887 ELSE IF ( ( t1 .GE. t_ref ) .AND. ( t1 .LT. -47. ) ) then !liq phas poor ES
3888 ew = es0 * exp(17.67 * t1 / ( t1 + 243.5))
3890 ELSE
3891 tk = t(i,k,j)
3892 rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + &
3893 c3 * (1. - tk / tf) + alog10(eis)
3894 ew = 10. ** rhs
3896 END IF
3898 ! Now sat vap pres obtained compute local vapor pressure
3900 ew = MAX ( ew , 0. ) * rh(i,k,j) * 0.01
3902 ! Now compute the specific humidity using the partial vapor
3903 ! pressures of water vapor (ew) and dry air (p-ew). The
3904 ! constants assume that the pressure is in hPa, so we divide
3905 ! the pressures by 100.
3907 q1 = mw_vap * ew
3908 q1 = q1 / (q1 + mw_air * (p(i,k,j)/100. - ew))
3910 q(i,k,j) = q1 / (1. - q1 )
3912 END IF
3914 END DO
3915 END DO
3916 END DO
3918 END IF
3920 END SUBROUTINE rh_to_mxrat
3922 !---------------------------------------------------------------------
3924 SUBROUTINE compute_eta ( znw , &
3925 eta_levels , max_eta , max_dz , &
3926 p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso , &
3927 ids , ide , jds , jde , kds , kde , &
3928 ims , ime , jms , jme , kms , kme , &
3929 its , ite , jts , jte , kts , kte )
3931 ! Compute eta levels, either using given values from the namelist (hardly
3932 ! a computation, yep, I know), or assuming a constant dz above the PBL,
3933 ! knowing p_top and the number of eta levels.
3935 IMPLICIT NONE
3937 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3938 ims , ime , jms , jme , kms , kme , &
3939 its , ite , jts , jte , kts , kte
3940 REAL , INTENT(IN) :: max_dz
3941 REAL , INTENT(IN) :: p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso
3942 INTEGER , INTENT(IN) :: max_eta
3943 REAL , DIMENSION (max_eta) , INTENT(IN) :: eta_levels
3945 REAL , DIMENSION (kts:kte) , INTENT(OUT) :: znw
3947 ! Local vars
3949 INTEGER :: k
3950 REAL :: mub , t_init , p_surf , pb, ztop, ztop_pbl , dz , temp
3951 REAL , DIMENSION(kts:kte) :: dnw
3953 INTEGER , PARAMETER :: prac_levels = 17
3954 INTEGER :: loop , loop1
3955 REAL , DIMENSION(prac_levels) :: znw_prac , znu_prac , dnw_prac
3956 REAL , DIMENSION(kts:kte) :: alb , phb
3958 ! Gee, do the eta levels come in from the namelist?
3960 IF ( ABS(eta_levels(1)+1.) .GT. 0.0000001 ) THEN
3962 ! Check to see if the array is oriented OK, we can easily fix an upside down oops.
3964 IF ( ( ABS(eta_levels(1 )-1.) .LT. 0.0000001 ) .AND. &
3965 ( ABS(eta_levels(kde)-0.) .LT. 0.0000001 ) ) THEN
3966 DO k = kds+1 , kde-1
3967 znw(k) = eta_levels(k)
3968 END DO
3969 znw( 1) = 1.
3970 znw(kde) = 0.
3971 ELSE IF ( ( ABS(eta_levels(kde)-1.) .LT. 0.0000001 ) .AND. &
3972 ( ABS(eta_levels(1 )-0.) .LT. 0.0000001 ) ) THEN
3973 DO k = kds+1 , kde-1
3974 znw(k) = eta_levels(kde+1-k)
3975 END DO
3976 znw( 1) = 1.
3977 znw(kde) = 0.
3978 ELSE
3979 CALL wrf_error_fatal ( 'First eta level should be 1.0 and the last 0.0 in namelist' )
3980 END IF
3982 ! Check to see if the input full-level eta array is monotonic.
3984 DO k = kds , kde-1
3985 IF ( znw(k) .LE. znw(k+1) ) THEN
3986 PRINT *,'eta on full levels is not monotonic'
3987 PRINT *,'eta (',k,') = ',znw(k)
3988 PRINT *,'eta (',k+1,') = ',znw(k+1)
3989 CALL wrf_error_fatal ( 'Fix non-monotonic "eta_levels" in the namelist.input file' )
3990 END IF
3991 END DO
3993 ! Compute eta levels assuming a constant delta z above the PBL.
3995 ELSE
3997 ! Compute top of the atmosphere with some silly levels. We just want to
3998 ! integrate to get a reasonable value for ztop. We use the planned PBL-esque
3999 ! levels, and then just coarse resolution above that. We know p_top, and we
4000 ! have the base state vars.
4002 p_surf = p00
4004 znw_prac = (/ 1.000 , 0.993 , 0.983 , 0.970 , 0.954 , 0.934 , 0.909 , &
4005 0.88 , 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /)
4007 DO k = 1 , prac_levels - 1
4008 znu_prac(k) = ( znw_prac(k) + znw_prac(k+1) ) * 0.5
4009 dnw_prac(k) = znw_prac(k+1) - znw_prac(k)
4010 END DO
4012 DO k = 1, prac_levels-1
4013 pb = znu_prac(k)*(p_surf - p_top) + p_top
4014 temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
4015 ! temp = t00 + A*LOG(pb/p00)
4016 t_init = temp*(p00/pb)**(r_d/cp) - t0
4017 alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
4018 END DO
4020 ! Base state mu is defined as base state surface pressure minus p_top
4022 mub = p_surf - p_top
4024 ! Integrate base geopotential, starting at terrain elevation.
4026 phb(1) = 0.
4027 DO k = 2,prac_levels
4028 phb(k) = phb(k-1) - dnw_prac(k-1)*mub*alb(k-1)
4029 END DO
4031 ! So, now we know the model top in meters. Get the average depth above the PBL
4032 ! of each of the remaining levels. We are going for a constant delta z thickness.
4034 ztop = phb(prac_levels) / g
4035 ztop_pbl = phb(8 ) / g
4036 dz = ( ztop - ztop_pbl ) / REAL ( kde - 8 )
4038 ! Standard levels near the surface so no one gets in trouble.
4040 DO k = 1 , 8
4041 znw(k) = znw_prac(k)
4042 END DO
4044 ! Using d phb(k)/ d eta(k) = -mub * alb(k), eqn 2.9
4045 ! Skamarock et al, NCAR TN 468. Use full levels, so
4046 ! use twice the thickness.
4048 DO k = 8, kte-1-2
4049 pb = znw(k) * (p_surf - p_top) + p_top
4050 temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
4051 ! temp = t00 + A*LOG(pb/p00)
4052 t_init = temp*(p00/pb)**(r_d/cp) - t0
4053 alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
4054 znw(k+1) = znw(k) - dz*g / ( mub*alb(k) )
4055 END DO
4056 znw(kte-2) = 0.000
4058 ! There is some iteration. We want the top level, ztop, to be
4059 ! consistent with the delta z, and we want the half level values
4060 ! to be consistent with the eta levels. The inner loop to 10 gets
4061 ! the eta levels very accurately, but has a residual at the top, due
4062 ! to dz changing. We reset dz five times, and then things seem OK.
4064 DO loop1 = 1 , 5
4065 DO loop = 1 , 10
4066 DO k = 8, kte-1-2
4067 pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
4068 temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
4069 ! temp = t00 + A*LOG(pb/p00)
4070 t_init = temp*(p00/pb)**(r_d/cp) - t0
4071 alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
4072 znw(k+1) = znw(k) - dz*g / ( mub*alb(k) )
4073 END DO
4074 IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN
4075 print *,'Converged znw(kte) should be about 0.0 = ',znw(kte-2)
4076 END IF
4077 znw(kte-2) = 0.000
4078 END DO
4080 ! Here is where we check the eta levels values we just computed.
4082 DO k = 1, kde-1-2
4083 pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
4084 temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
4085 ! temp = t00 + A*LOG(pb/p00)
4086 t_init = temp*(p00/pb)**(r_d/cp) - t0
4087 alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
4088 END DO
4090 phb(1) = 0.
4091 DO k = 2,kde-2
4092 phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1)
4093 END DO
4095 ! Reset the model top and the dz, and iterate.
4097 ztop = phb(kde-2)/g
4098 ztop_pbl = phb(8)/g
4099 dz = ( ztop - ztop_pbl ) / REAL ( (kde-2) - 8 )
4100 END DO
4102 IF ( dz .GT. max_dz ) THEN
4103 print *,'z (m) = ',phb(1)/g
4104 do k = 2 ,kte-2
4105 print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g
4106 end do
4107 print *,'dz (m) above fixed eta levels = ',dz
4108 print *,'namelist max_dz (m) = ',max_dz
4109 print *,'namelist p_top (Pa) = ',p_top
4110 CALL wrf_debug ( 0, 'You need one of three things:' )
4111 CALL wrf_debug ( 0, '1) More eta levels to reduce the dz: e_vert' )
4112 CALL wrf_debug ( 0, '2) A lower p_top so your total height is reduced: p_top_requested')
4113 CALL wrf_debug ( 0, '3) Increase the maximum allowable eta thickness: max_dz')
4114 CALL wrf_debug ( 0, 'All are namelist options')
4115 CALL wrf_error_fatal ( 'dz above fixed eta levels is too large')
4116 END IF
4118 ! Add those 2 levels back into the middle, just above the 8 levels
4119 ! that semi define a boundary layer. After we open up the levels,
4120 ! then we just linearly interpolate in znw. So now levels 1-8 are
4121 ! specified as the fixed boundary layer levels given in this routine.
4122 ! The top levels, 12 through kte are those computed. The middle
4123 ! levels 9, 10, and 11 are equi-spaced in znw, and are each 1/2 the
4124 ! the znw thickness of levels 11 through 12.
4126 DO k = kte-2 , 9 , -1
4127 znw(k+2) = znw(k)
4128 END DO
4130 znw( 9) = 0.75 * znw( 8) + 0.25 * znw(12)
4131 znw(10) = 0.50 * znw( 8) + 0.50 * znw(12)
4132 znw(11) = 0.25 * znw( 8) + 0.75 * znw(12)
4134 END IF
4136 END SUBROUTINE compute_eta
4138 !---------------------------------------------------------------------
4140 SUBROUTINE monthly_min_max ( field_in , field_min , field_max , &
4141 ids , ide , jds , jde , kds , kde , &
4142 ims , ime , jms , jme , kms , kme , &
4143 its , ite , jts , jte , kts , kte )
4145 ! Plow through each month, find the max, min values for each i,j.
4147 IMPLICIT NONE
4149 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4150 ims , ime , jms , jme , kms , kme , &
4151 its , ite , jts , jte , kts , kte
4153 REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in
4154 REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max
4156 ! Local vars
4158 INTEGER :: i , j , l
4159 REAL :: minner , maxxer
4161 DO j = jts , MIN(jde-1,jte)
4162 DO i = its , MIN(ide-1,ite)
4163 minner = field_in(i,1,j)
4164 maxxer = field_in(i,1,j)
4165 DO l = 2 , 12
4166 IF ( field_in(i,l,j) .LT. minner ) THEN
4167 minner = field_in(i,l,j)
4168 END IF
4169 IF ( field_in(i,l,j) .GT. maxxer ) THEN
4170 maxxer = field_in(i,l,j)
4171 END IF
4172 END DO
4173 field_min(i,j) = minner
4174 field_max(i,j) = maxxer
4175 END DO
4176 END DO
4178 END SUBROUTINE monthly_min_max
4180 !---------------------------------------------------------------------
4182 SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , &
4183 ids , ide , jds , jde , kds , kde , &
4184 ims , ime , jms , jme , kms , kme , &
4185 its , ite , jts , jte , kts , kte )
4187 ! Linrarly in time interpolate data to a current valid time. The data is
4188 ! assumed to come in "monthly", valid at the 15th of every month.
4190 IMPLICIT NONE
4192 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4193 ims , ime , jms , jme , kms , kme , &
4194 its , ite , jts , jte , kts , kte
4196 CHARACTER (LEN=24) , INTENT(IN) :: date_str
4197 REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in
4198 REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out
4200 ! Local vars
4202 INTEGER :: i , j , l
4203 INTEGER , DIMENSION(0:13) :: middle
4204 INTEGER :: target_julyr , target_julday , target_date
4205 INTEGER :: julyr , julday , int_month , month1 , month2
4206 REAL :: gmt
4207 CHARACTER (LEN=4) :: yr
4208 CHARACTER (LEN=2) :: mon , day15
4211 WRITE(day15,FMT='(I2.2)') 15
4212 DO l = 1 , 12
4213 WRITE(mon,FMT='(I2.2)') l
4214 CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt )
4215 middle(l) = julyr*1000 + julday
4216 END DO
4218 l = 0
4219 middle(l) = middle( 1) - 31
4221 l = 13
4222 middle(l) = middle(12) + 31
4224 CALL get_julgmt ( date_str , target_julyr , target_julday , gmt )
4225 target_date = target_julyr * 1000 + target_julday
4226 find_month : DO l = 0 , 12
4227 IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN
4228 DO j = jts , MIN ( jde-1 , jte )
4229 DO i = its , MIN (ide-1 , ite )
4230 int_month = l
4231 IF ( ( int_month .EQ. 0 ) .OR. ( int_month .EQ. 12 ) ) THEN
4232 month1 = 12
4233 month2 = 1
4234 ELSE
4235 month1 = int_month
4236 month2 = month1 + 1
4237 END IF
4238 field_out(i,j) = ( field_in(i,month2,j) * ( target_date - middle(l) ) + &
4239 field_in(i,month1,j) * ( middle(l+1) - target_date ) ) / &
4240 ( middle(l+1) - middle(l) )
4241 END DO
4242 END DO
4243 EXIT find_month
4244 END IF
4245 END DO find_month
4247 END SUBROUTINE monthly_interp_to_date
4249 !---------------------------------------------------------------------
4251 SUBROUTINE sfcprs (t, q, height, pslv, ter, avgsfct, p, &
4252 psfc, ez_method, &
4253 ids , ide , jds , jde , kds , kde , &
4254 ims , ime , jms , jme , kms , kme , &
4255 its , ite , jts , jte , kts , kte )
4258 ! Computes the surface pressure using the input height,
4259 ! temperature and q (already computed from relative
4260 ! humidity) on p surfaces. Sea level pressure is used
4261 ! to extrapolate a first guess.
4263 IMPLICIT NONE
4265 REAL, PARAMETER :: g = 9.8
4266 REAL, PARAMETER :: gamma = 6.5E-3
4267 REAL, PARAMETER :: pconst = 10000.0
4268 REAL, PARAMETER :: Rd = 287.
4269 REAL, PARAMETER :: TC = 273.15 + 17.5
4271 REAL, PARAMETER :: gammarg = gamma * Rd / g
4272 REAL, PARAMETER :: rov2 = Rd / 2.
4274 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4275 ims , ime , jms , jme , kms , kme , &
4276 its , ite , jts , jte , kts , kte
4277 LOGICAL , INTENT ( IN ) :: ez_method
4279 REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p
4280 REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: pslv , ter, avgsfct
4281 REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
4283 INTEGER :: i
4284 INTEGER :: j
4285 INTEGER :: k
4286 INTEGER , DIMENSION (its:ite,jts:jte) :: k500 , k700 , k850
4288 LOGICAL :: l1
4289 LOGICAL :: l2
4290 LOGICAL :: l3
4291 LOGICAL :: OK
4293 REAL :: gamma78 ( its:ite,jts:jte )
4294 REAL :: gamma57 ( its:ite,jts:jte )
4295 REAL :: ht ( its:ite,jts:jte )
4296 REAL :: p1 ( its:ite,jts:jte )
4297 REAL :: t1 ( its:ite,jts:jte )
4298 REAL :: t500 ( its:ite,jts:jte )
4299 REAL :: t700 ( its:ite,jts:jte )
4300 REAL :: t850 ( its:ite,jts:jte )
4301 REAL :: tfixed ( its:ite,jts:jte )
4302 REAL :: tsfc ( its:ite,jts:jte )
4303 REAL :: tslv ( its:ite,jts:jte )
4305 ! We either compute the surface pressure from a time averaged surface temperature
4306 ! (what we will call the "easy way"), or we try to remove the diurnal impact on the
4307 ! surface temperature (what we will call the "other way"). Both are essentially
4308 ! corrections to a sea level pressure with a high-resolution topography field.
4310 IF ( ez_method ) THEN
4312 DO j = jts , MIN(jde-1,jte)
4313 DO i = its , MIN(ide-1,ite)
4314 psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / avgsfct(i,j) ) ** ( - g / ( Rd * gamma ) )
4315 END DO
4316 END DO
4318 ELSE
4320 ! Find the locations of the 850, 700 and 500 mb levels.
4322 k850 = 0 ! find k at: P=850
4323 k700 = 0 ! P=700
4324 k500 = 0 ! P=500
4326 i = its
4327 j = jts
4328 DO k = kts+1 , kte
4329 IF (NINT(p(i,k,j)) .EQ. 85000) THEN
4330 k850(i,j) = k
4331 ELSE IF (NINT(p(i,k,j)) .EQ. 70000) THEN
4332 k700(i,j) = k
4333 ELSE IF (NINT(p(i,k,j)) .EQ. 50000) THEN
4334 k500(i,j) = k
4335 END IF
4336 END DO
4338 IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN
4340 DO j = jts , MIN(jde-1,jte)
4341 DO i = its , MIN(ide-1,ite)
4342 psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / t(i,1,j) ) ** ( - g / ( Rd * gamma ) )
4343 END DO
4344 END DO
4346 RETURN
4347 #if 0
4349 ! Possibly it is just that we have a generalized vertical coord, so we do not
4350 ! have the values exactly. Do a simple assignment to a close vertical level.
4352 DO j = jts , MIN(jde-1,jte)
4353 DO i = its , MIN(ide-1,ite)
4354 DO k = kts+1 , kte-1
4355 IF ( ( p(i,k,j) - 85000. ) * ( p(i,k+1,j) - 85000. ) .LE. 0.0 ) THEN
4356 k850(i,j) = k
4357 END IF
4358 IF ( ( p(i,k,j) - 70000. ) * ( p(i,k+1,j) - 70000. ) .LE. 0.0 ) THEN
4359 k700(i,j) = k
4360 END IF
4361 IF ( ( p(i,k,j) - 50000. ) * ( p(i,k+1,j) - 50000. ) .LE. 0.0 ) THEN
4362 k500(i,j) = k
4363 END IF
4364 END DO
4365 END DO
4366 END DO
4368 ! If we *still* do not have the k levels, punt. I mean, we did try.
4370 OK = .TRUE.
4371 DO j = jts , MIN(jde-1,jte)
4372 DO i = its , MIN(ide-1,ite)
4373 IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN
4374 OK = .FALSE.
4375 PRINT '(A)','(i,j) = ',i,j,' Error in finding p level for 850, 700 or 500 hPa.'
4376 DO K = kts+1 , kte
4377 PRINT '(A,I3,A,F10.2,A)','K = ',k,' PRESSURE = ',p(i,k,j),' Pa'
4378 END DO
4379 PRINT '(A)','Expected 850, 700, and 500 mb values, at least.'
4380 END IF
4381 END DO
4382 END DO
4383 IF ( .NOT. OK ) THEN
4384 CALL wrf_error_fatal ( 'wrong pressure levels' )
4385 END IF
4386 #endif
4388 ! We are here if the data is isobaric and we found the levels for 850, 700,
4389 ! and 500 mb right off the bat.
4391 ELSE
4392 DO j = jts , MIN(jde-1,jte)
4393 DO i = its , MIN(ide-1,ite)
4394 k850(i,j) = k850(its,jts)
4395 k700(i,j) = k700(its,jts)
4396 k500(i,j) = k500(its,jts)
4397 END DO
4398 END DO
4399 END IF
4401 ! The 850 hPa level of geopotential height is called something special.
4403 DO j = jts , MIN(jde-1,jte)
4404 DO i = its , MIN(ide-1,ite)
4405 ht(i,j) = height(i,k850(i,j),j)
4406 END DO
4407 END DO
4409 ! The variable ht is now -ter/ht(850 hPa). The plot thickens.
4411 DO j = jts , MIN(jde-1,jte)
4412 DO i = its , MIN(ide-1,ite)
4413 ht(i,j) = -ter(i,j) / ht(i,j)
4414 END DO
4415 END DO
4417 ! Make an isothermal assumption to get a first guess at the surface
4418 ! pressure. This is to tell us which levels to use for the lapse
4419 ! rates in a bit.
4421 DO j = jts , MIN(jde-1,jte)
4422 DO i = its , MIN(ide-1,ite)
4423 psfc(i,j) = pslv(i,j) * (pslv(i,j) / p(i,k850(i,j),j)) ** ht(i,j)
4424 END DO
4425 END DO
4427 ! Get a pressure more than pconst Pa above the surface - p1. The
4428 ! p1 is the top of the level that we will use for our lapse rate
4429 ! computations.
4431 DO j = jts , MIN(jde-1,jte)
4432 DO i = its , MIN(ide-1,ite)
4433 IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN
4434 p1(i,j) = 85000.
4435 ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0. ) THEN
4436 p1(i,j) = psfc(i,j) - pconst
4437 ELSE
4438 p1(i,j) = 50000.
4439 END IF
4440 END DO
4441 END DO
4443 ! Compute virtual temperatures for k850, k700, and k500 layers. Now
4444 ! you see why we wanted Q on pressure levels, it all is beginning
4445 ! to make sense.
4447 DO j = jts , MIN(jde-1,jte)
4448 DO i = its , MIN(ide-1,ite)
4449 t850(i,j) = t(i,k850(i,j),j) * (1. + 0.608 * q(i,k850(i,j),j))
4450 t700(i,j) = t(i,k700(i,j),j) * (1. + 0.608 * q(i,k700(i,j),j))
4451 t500(i,j) = t(i,k500(i,j),j) * (1. + 0.608 * q(i,k500(i,j),j))
4452 END DO
4453 END DO
4455 ! Compute lapse rates between these three levels. These are
4456 ! environmental values for each (i,j).
4458 DO j = jts , MIN(jde-1,jte)
4459 DO i = its , MIN(ide-1,ite)
4460 gamma78(i,j) = ALOG(t850(i,j) / t700(i,j)) / ALOG (p(i,k850(i,j),j) / p(i,k700(i,j),j) )
4461 gamma57(i,j) = ALOG(t700(i,j) / t500(i,j)) / ALOG (p(i,k700(i,j),j) / p(i,k500(i,j),j) )
4462 END DO
4463 END DO
4465 DO j = jts , MIN(jde-1,jte)
4466 DO i = its , MIN(ide-1,ite)
4467 IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN
4468 t1(i,j) = t850(i,j)
4469 ELSE IF ( ( psfc(i,j) - 85000. ) .GE. 0. ) THEN
4470 t1(i,j) = t700(i,j) * (p1(i,j) / (p(i,k700(i,j),j))) ** gamma78(i,j)
4471 ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0.) THEN
4472 t1(i,j) = t500(i,j) * (p1(i,j) / (p(i,k500(i,j),j))) ** gamma57(i,j)
4473 ELSE
4474 t1(i,j) = t500(i,j)
4475 ENDIF
4476 END DO
4477 END DO
4479 ! From our temperature way up in the air, we extrapolate down to
4480 ! the sea level to get a guess at the sea level temperature.
4482 DO j = jts , MIN(jde-1,jte)
4483 DO i = its , MIN(ide-1,ite)
4484 tslv(i,j) = t1(i,j) * (pslv(i,j) / p1(i,j)) ** gammarg
4485 END DO
4486 END DO
4488 ! The new surface temperature is computed from the with new sea level
4489 ! temperature, just using the elevation and a lapse rate. This lapse
4490 ! rate is -6.5 K/km.
4492 DO j = jts , MIN(jde-1,jte)
4493 DO i = its , MIN(ide-1,ite)
4494 tsfc(i,j) = tslv(i,j) - gamma * ter(i,j)
4495 END DO
4496 END DO
4498 ! A small correction to the sea-level temperature, in case it is too warm.
4500 DO j = jts , MIN(jde-1,jte)
4501 DO i = its , MIN(ide-1,ite)
4502 tfixed(i,j) = tc - 0.005 * (tsfc(i,j) - tc) ** 2
4503 END DO
4504 END DO
4506 DO j = jts , MIN(jde-1,jte)
4507 DO i = its , MIN(ide-1,ite)
4508 l1 = tslv(i,j) .LT. tc
4509 l2 = tsfc(i,j) .LE. tc
4510 l3 = .NOT. l1
4511 IF ( l2 .AND. l3 ) THEN
4512 tslv(i,j) = tc
4513 ELSE IF ( ( .NOT. l2 ) .AND. l3 ) THEN
4514 tslv(i,j) = tfixed(i,j)
4515 END IF
4516 END DO
4517 END DO
4519 ! Finally, we can get to the surface pressure.
4521 DO j = jts , MIN(jde-1,jte)
4522 DO i = its , MIN(ide-1,ite)
4523 p1(i,j) = - ter(i,j) * g / ( rov2 * ( tsfc(i,j) + tslv(i,j) ) )
4524 psfc(i,j) = pslv(i,j) * EXP ( p1(i,j) )
4525 END DO
4526 END DO
4528 END IF
4530 ! Surface pressure and sea-level pressure are the same at sea level.
4532 ! DO j = jts , MIN(jde-1,jte)
4533 ! DO i = its , MIN(ide-1,ite)
4534 ! IF ( ABS ( ter(i,j) ) .LT. 0.1 ) THEN
4535 ! psfc(i,j) = pslv(i,j)
4536 ! END IF
4537 ! END DO
4538 ! END DO
4540 END SUBROUTINE sfcprs
4542 !---------------------------------------------------------------------
4544 SUBROUTINE sfcprs2(t, q, height, psfc_in, ter, avgsfct, p, &
4545 psfc, ez_method, &
4546 ids , ide , jds , jde , kds , kde , &
4547 ims , ime , jms , jme , kms , kme , &
4548 its , ite , jts , jte , kts , kte )
4551 ! Computes the surface pressure using the input height,
4552 ! temperature and q (already computed from relative
4553 ! humidity) on p surfaces. Sea level pressure is used
4554 ! to extrapolate a first guess.
4556 IMPLICIT NONE
4558 REAL, PARAMETER :: g = 9.8
4559 REAL, PARAMETER :: Rd = 287.
4561 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4562 ims , ime , jms , jme , kms , kme , &
4563 its , ite , jts , jte , kts , kte
4564 LOGICAL , INTENT ( IN ) :: ez_method
4566 REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p
4567 REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: psfc_in , ter, avgsfct
4568 REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
4570 INTEGER :: i
4571 INTEGER :: j
4572 INTEGER :: k
4574 REAL :: tv_sfc_avg , tv_sfc , del_z
4576 ! Compute the new surface pressure from the old surface pressure, and a
4577 ! known change in elevation at the surface.
4579 ! del_z = diff in surface topo, lo-res vs hi-res
4580 ! psfc = psfc_in * exp ( g del_z / (Rd Tv_sfc ) )
4583 IF ( ez_method ) THEN
4584 DO j = jts , MIN(jde-1,jte)
4585 DO i = its , MIN(ide-1,ite)
4586 tv_sfc_avg = avgsfct(i,j) * (1. + 0.608 * q(i,1,j))
4587 del_z = height(i,1,j) - ter(i,j)
4588 psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc_avg ) )
4589 END DO
4590 END DO
4591 ELSE
4592 DO j = jts , MIN(jde-1,jte)
4593 DO i = its , MIN(ide-1,ite)
4594 tv_sfc = t(i,1,j) * (1. + 0.608 * q(i,1,j))
4595 del_z = height(i,1,j) - ter(i,j)
4596 psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc ) )
4597 END DO
4598 END DO
4599 END IF
4601 END SUBROUTINE sfcprs2
4603 !---------------------------------------------------------------------
4605 SUBROUTINE sfcprs3( height , p , ter , slp , psfc , &
4606 ids , ide , jds , jde , kds , kde , &
4607 ims , ime , jms , jme , kms , kme , &
4608 its , ite , jts , jte , kts , kte )
4610 ! Computes the surface pressure by vertically interpolating
4611 ! linearly (or log) in z the pressure, to the targeted topography.
4613 IMPLICIT NONE
4615 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4616 ims , ime , jms , jme , kms , kme , &
4617 its , ite , jts , jte , kts , kte
4619 REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: height, p
4620 REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: ter , slp
4621 REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
4623 INTEGER :: i
4624 INTEGER :: j
4625 INTEGER :: k
4627 LOGICAL :: found_loc
4629 REAL :: zl , zu , pl , pu , zm
4631 ! Loop over each grid point
4633 DO j = jts , MIN(jde-1,jte)
4634 DO i = its , MIN(ide-1,ite)
4636 ! Find the trapping levels
4638 found_loc = .FALSE.
4640 ! Normal sort of scenario - the model topography is somewhere between
4641 ! the height values of 1000 mb and the top of the model.
4643 found_k_loc : DO k = kts+1 , kte-2
4644 IF ( ( height(i,k ,j) .LE. ter(i,j) ) .AND. &
4645 ( height(i,k+1,j) .GT. ter(i,j) ) ) THEN
4646 zl = height(i,k ,j)
4647 zu = height(i,k+1,j)
4648 zm = ter(i,j)
4649 pl = p(i,k ,j)
4650 pu = p(i,k+1,j)
4651 psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
4652 found_loc = .TRUE.
4653 EXIT found_k_loc
4654 END IF
4655 END DO found_k_loc
4657 ! Interpolate betwixt slp and the first isobaric level above - this is probably the
4658 ! usual thing over the ocean.
4660 IF ( .NOT. found_loc ) THEN
4661 IF ( slp(i,j) .GE. p(i,2,j) ) THEN
4662 zl = 0.
4663 zu = height(i,2 ,j)
4664 zm = ter(i,j)
4665 pl = slp(i,j)
4666 pu = p(i,2 ,j)
4667 psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
4668 found_loc = .TRUE.
4669 ELSE
4670 found_slp_loc : DO k = kts+1 , kte-2
4671 IF ( ( slp(i,j) .GE. p(i,k+1,j) ) .AND. &
4672 ( slp(i,j) .LT. p(i,k ,j) ) ) THEN
4673 zl = 0.
4674 zu = height(i,k+1,j)
4675 zm = ter(i,j)
4676 pl = slp(i,j)
4677 pu = p(i,k+1,j)
4678 psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
4679 found_loc = .TRUE.
4680 EXIT found_slp_loc
4681 END IF
4682 END DO found_slp_loc
4683 END IF
4684 END IF
4686 ! Did we do what we wanted done.
4688 IF ( .NOT. found_loc ) THEN
4689 print *,'i,j = ',i,j
4690 print *,'p column = ',p(i,2:,j)
4691 print *,'z column = ',height(i,2:,j)
4692 print *,'model topo = ',ter(i,j)
4693 CALL wrf_error_fatal ( ' probs with sfc p computation ' )
4694 END IF
4696 END DO
4697 END DO
4699 END SUBROUTINE sfcprs3
4700 !---------------------------------------------------------------------
4702 SUBROUTINE filter_topo ( ht_in , xlat , msftx , fft_filter_lat , &
4703 ids , ide , jds , jde , kds , kde , &
4704 ims , ime , jms , jme , kms , kme , &
4705 its , ite , jts , jte , kts , kte )
4707 IMPLICIT NONE
4709 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4710 ims , ime , jms , jme , kms , kme , &
4711 its , ite , jts , jte , kts , kte
4713 REAL , INTENT(IN) :: fft_filter_lat
4714 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: ht_in
4715 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: xlat , msftx
4718 ! Local vars
4720 INTEGER :: i , j , j_lat_pos , j_lat_neg
4721 INTEGER :: i_kicker , ik , i1, i2, i3, i4
4722 REAL :: length_scale , sum
4723 REAL , DIMENSION(its:ite,jts:jte) :: ht_out
4725 ! The filtering is a simple average on a latitude loop. Possibly a LONG list of
4726 ! numbers. We assume that ALL of the 2d arrays have been transposed so that
4727 ! each patch has the entire domain size of the i-dim local.
4729 IF ( ( its .NE. ids ) .OR. ( ite .NE. ide ) ) THEN
4730 CALL wrf_error_fatal ( 'filtering assumes all values on X' )
4731 END IF
4733 ! Starting at the south pole, we find where the
4734 ! grid distance is big enough, then go back a point. Continuing to the
4735 ! north pole, we find the first small grid distance. These are the
4736 ! computational latitude loops and the associated computational poles.
4738 j_lat_neg = 0
4739 j_lat_pos = jde + 1
4740 loop_neg : DO j = jts , MIN(jde-1,jte)
4741 IF ( xlat(its,j) .LT. 0.0 ) THEN
4742 IF ( ABS(xlat(its,j)) .LT. fft_filter_lat ) THEN
4743 j_lat_neg = j - 1
4744 EXIT loop_neg
4745 END IF
4746 END IF
4747 END DO loop_neg
4749 loop_pos : DO j = jts , MIN(jde-1,jte)
4750 IF ( xlat(its,j) .GT. 0.0 ) THEN
4751 IF ( xlat(its,j) .GE. fft_filter_lat ) THEN
4752 j_lat_pos = j
4753 EXIT loop_pos
4754 END IF
4755 END IF
4756 END DO loop_pos
4758 ! Set output values to initial input topo values for whole patch.
4760 DO j = jts , MIN(jde-1,jte)
4761 DO i = its , MIN(ide-1,ite)
4762 ht_out(i,j) = ht_in(i,j)
4763 END DO
4764 END DO
4766 ! Filter the topo at the negative lats.
4768 DO j = j_lat_neg , jts , -1
4769 i_kicker = MIN( MAX ( NINT(msftx(its,j)) , 1 ) , (ide - ids) / 2 )
4770 print *,'j = ' , j, ', kicker = ',i_kicker
4771 DO i = its , MIN(ide-1,ite)
4772 IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
4773 sum = 0.0
4774 DO ik = 1 , i_kicker
4775 sum = sum + ht_in(i+ik,j) + ht_in(i-ik,j)
4776 END DO
4777 ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
4778 ELSE IF ( ( i - i_kicker .LT. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
4779 sum = 0.0
4780 DO ik = 1 , i_kicker
4781 sum = sum + ht_in(i+ik,j)
4782 END DO
4783 i1 = i - i_kicker + ide -1
4784 i2 = ide-1
4785 i3 = ids
4786 i4 = i-1
4787 DO ik = i1 , i2
4788 sum = sum + ht_in(ik,j)
4789 END DO
4790 DO ik = i3 , i4
4791 sum = sum + ht_in(ik,j)
4792 END DO
4793 ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
4794 ELSE IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .GT. ide-1 ) ) THEN
4795 sum = 0.0
4796 DO ik = 1 , i_kicker
4797 sum = sum + ht_in(i-ik,j)
4798 END DO
4799 i1 = i+1
4800 i2 = ide-1
4801 i3 = ids
4802 i4 = ids + ( i_kicker+i ) - ide
4803 DO ik = i1 , i2
4804 sum = sum + ht_in(ik,j)
4805 END DO
4806 DO ik = i3 , i4
4807 sum = sum + ht_in(ik,j)
4808 END DO
4809 ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
4810 END IF
4811 END DO
4812 END DO
4814 ! Filter the topo at the positive lats.
4816 DO j = j_lat_pos , MIN(jde-1,jte)
4817 i_kicker = MIN( MAX ( NINT(msftx(its,j)) , 1 ) , (ide - ids) / 2 )
4818 print *,'j = ' , j, ', kicker = ',i_kicker
4819 DO i = its , MIN(ide-1,ite)
4820 IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
4821 sum = 0.0
4822 DO ik = 1 , i_kicker
4823 sum = sum + ht_in(i+ik,j) + ht_in(i-ik,j)
4824 END DO
4825 ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
4826 ELSE IF ( ( i - i_kicker .LT. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
4827 sum = 0.0
4828 DO ik = 1 , i_kicker
4829 sum = sum + ht_in(i+ik,j)
4830 END DO
4831 i1 = i - i_kicker + ide -1
4832 i2 = ide-1
4833 i3 = ids
4834 i4 = i-1
4835 DO ik = i1 , i2
4836 sum = sum + ht_in(ik,j)
4837 END DO
4838 DO ik = i3 , i4
4839 sum = sum + ht_in(ik,j)
4840 END DO
4841 ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
4842 ELSE IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .GT. ide-1 ) ) THEN
4843 sum = 0.0
4844 DO ik = 1 , i_kicker
4845 sum = sum + ht_in(i-ik,j)
4846 END DO
4847 i1 = i+1
4848 i2 = ide-1
4849 i3 = ids
4850 i4 = ids + ( i_kicker+i ) - ide
4851 DO ik = i1 , i2
4852 sum = sum + ht_in(ik,j)
4853 END DO
4854 DO ik = i3 , i4
4855 sum = sum + ht_in(ik,j)
4856 END DO
4857 ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
4858 END IF
4859 END DO
4860 END DO
4862 ! Set output values to initial input topo values for whole patch.
4864 DO j = jts , MIN(jde-1,jte)
4865 DO i = its , MIN(ide-1,ite)
4866 ht_in(i,j) = ht_out(i,j)
4867 END DO
4868 END DO
4870 END SUBROUTINE filter_topo
4872 !---------------------------------------------------------------------
4874 SUBROUTINE init_module_initialize
4875 END SUBROUTINE init_module_initialize
4877 !---------------------------------------------------------------------
4879 END MODULE module_initialize_real
4880 #endif