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