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