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