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