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