Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

wrf_timeseries.F @ 1419

Last change on this file since 1419 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 19.9 KB

Line
1	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2	! This routine prints out the current value of variables at all specified
3	! time series locations that are within the current patch.
4	!
5	! Michael G. Duda -- 25 August 2005
6	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
7	SUBROUTINE calc_ts_locations( grid )
8
9	USE module_domain, ONLY : domain, get_ijk_from_grid
10	USE module_configure, ONLY : model_config_rec, grid_config_rec_type, model_to_grid_config_rec
11	USE module_dm, ONLY : wrf_dm_min_real
12	USE module_llxy
13	USE module_state_description
14
15	IMPLICIT NONE
16
17	! Arguments
18	TYPE (domain), INTENT(INOUT) :: grid
19
20	! Externals
21	LOGICAL, EXTERNAL :: wrf_dm_on_monitor
22	INTEGER, EXTERNAL :: get_unused_unit
23
24	! Local variables
25	INTEGER :: ntsloc_temp
26	INTEGER :: i, k, iunit
27	REAL :: ts_rx, ts_ry, ts_xlat, ts_xlong, ts_hgt
28	REAL :: known_lat, known_lon
29	CHARACTER (LEN=132) :: message
30	TYPE (PROJ_INFO) :: ts_proj
31	TYPE (grid_config_rec_type) :: config_flags
32
33	INTEGER :: ids, ide, jds, jde, kds, kde, &
34	ims, ime, jms, jme, kms, kme, &
35	ips, ipe, jps, jpe, kps, kpe, &
36	imsx, imex, jmsx, jmex, kmsx, kmex, &
37	ipsx, ipex, jpsx, jpex, kpsx, kpex, &
38	imsy, imey, jmsy, jmey, kmsy, kmey, &
39	ipsy, ipey, jpsy, jpey, kpsy, kpey
40
41
42	IF ( grid%ntsloc .LE. 0 ) RETURN
43
44	#if ((EM_CORE == 1) && (DA_CORE != 1))
45	IF ( grid%dfi_stage == DFI_FST ) THEN
46	#endif
47	CALL get_ijk_from_grid ( grid , &
48	ids, ide, jds, jde, kds, kde, &
49	ims, ime, jms, jme, kms, kme, &
50	ips, ipe, jps, jpe, kps, kpe, &
51	imsx, imex, jmsx, jmex, kmsx, kmex, &
52	ipsx, ipex, jpsx, jpex, kpsx, kpex, &
53	imsy, imey, jmsy, jmey, kmsy, kmey, &
54	ipsy, ipey, jpsy, jpey, kpsy, kpey )
55
56	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
57
58	! Set up map transformation structure
59	CALL map_init(ts_proj)
60
61	IF (ips <= 1 .AND. 1 <= ipe .AND. &
62	jps <= 1 .AND. 1 <= jpe) THEN
63	known_lat = grid%xlat(1,1)
64	known_lon = grid%xlong(1,1)
65	ELSE
66	known_lat = 9999.
67	known_lon = 9999.
68	END IF
69	known_lat = wrf_dm_min_real(known_lat)
70	known_lon = wrf_dm_min_real(known_lon)
71
72	! Mercator
73	IF (config_flags%map_proj == PROJ_MERC) THEN
74	CALL map_set(PROJ_MERC, ts_proj, &
75	truelat1 = config_flags%truelat1, &
76	lat1 = known_lat, &
77	lon1 = known_lon, &
78	knowni = 1., &
79	knownj = 1., &
80	dx = config_flags%dx)
81
82	! Lambert conformal
83	ELSE IF (config_flags%map_proj == PROJ_LC) THEN
84	CALL map_set(PROJ_LC, ts_proj, &
85	truelat1 = config_flags%truelat1, &
86	truelat2 = config_flags%truelat2, &
87	stdlon = config_flags%stand_lon, &
88	lat1 = known_lat, &
89	lon1 = known_lon, &
90	knowni = 1., &
91	knownj = 1., &
92	dx = config_flags%dx)
93
94	! Polar stereographic
95	ELSE IF (config_flags%map_proj == PROJ_PS) THEN
96	CALL map_set(PROJ_PS, ts_proj, &
97	truelat1 = config_flags%truelat1, &
98	stdlon = config_flags%stand_lon, &
99	lat1 = known_lat, &
100	lon1 = known_lon, &
101	knowni = 1., &
102	knownj = 1., &
103	dx = config_flags%dx)
104
105	#if (EM_CORE == 1)
106	! Cassini (global ARW)
107	ELSE IF (config_flags%map_proj == PROJ_CASSINI) THEN
108	CALL map_set(PROJ_CASSINI, ts_proj, &
109	latinc = grid%dy360.0/(2.0EARTH_RADIUS_M*PI), &
110	loninc = grid%dx360.0/(2.0EARTH_RADIUS_M*PI), &
111	lat1 = known_lat, &
112	lon1 = known_lon, &
113	lat0 = config_flags%pole_lat, &
114	lon0 = config_flags%pole_lon, &
115	knowni = 1., &
116	knownj = 1., &
117	stdlon = config_flags%stand_lon)
118	#endif
119
120	! Rotated latitude-longitude
121	ELSE IF (config_flags%map_proj == PROJ_ROTLL) THEN
122	CALL map_set(PROJ_ROTLL, ts_proj, &
123	! I have no idea how this should work for NMM nested domains
124	ixdim = grid%e_we-1, &
125	jydim = grid%e_sn-1, &
126	phi = real(grid%e_sn-2)*grid%dy/2.0, &
127	lambda = real(grid%e_we-2)*grid%dx, &
128	lat1 = config_flags%cen_lat, &
129	lon1 = config_flags%cen_lon, &
130	latinc = grid%dy, &
131	loninc = grid%dx, &
132	stagger = HH)
133
134	END IF
135
136	! Determine time series locations for domain
137	IF (.NOT. grid%have_calculated_tslocs) THEN
138	grid%have_calculated_tslocs = .TRUE.
139	WRITE(message, '(A43,I3)') 'Computing time series locations for domain ', grid%id
140	CALL wrf_message(message)
141
142	ntsloc_temp = 0
143	DO k=1,grid%ntsloc
144
145	IF (config_flags%map_proj == 0) THEN ! For idealized cases, no map transformation needed
146	ts_rx = grid%lattsloc(k) ! NB: (x,y) = (lat,lon) rather than (x,y) = (lon,lat)
147	ts_ry = grid%lontsloc(k)
148	ELSE
149	CALL latlon_to_ij(ts_proj, grid%lattsloc(k), grid%lontsloc(k), ts_rx, ts_ry)
150	END IF
151
152
153	ntsloc_temp = ntsloc_temp + 1
154	grid%itsloc(ntsloc_temp) = NINT(ts_rx)
155	grid%jtsloc(ntsloc_temp) = NINT(ts_ry)
156	grid%id_tsloc(ntsloc_temp) = k
157
158	! Is point outside of domain (or on the edge of domain)?
159	IF (grid%itsloc(ntsloc_temp) < ids .OR. grid%itsloc(ntsloc_temp) > ide .OR. &
160	grid%jtsloc(ntsloc_temp) < jds .OR. grid%jtsloc(ntsloc_temp) > jde) THEN
161	ntsloc_temp = ntsloc_temp - 1
162
163	END IF
164
165	END DO
166
167	grid%next_ts_time = 1
168
169	grid%ntsloc_domain = ntsloc_temp
170
171	DO k=1,grid%ntsloc_domain
172
173	! If location is outside of patch, we need to get lat/lon of TS grid cell from another patch
174	IF (grid%itsloc(k) < ips .OR. grid%itsloc(k) > ipe .OR. &
175	grid%jtsloc(k) < jps .OR. grid%jtsloc(k) > jpe) THEN
176	ts_xlat = 1.E30
177	ts_xlong = 1.E30
178	ts_hgt = 1.E30
179	ELSE
180	ts_xlat = grid%xlat(grid%itsloc(k),grid%jtsloc(k))
181	ts_xlong = grid%xlong(grid%itsloc(k),grid%jtsloc(k))
182	#if (EM_CORE == 1)
183	ts_hgt = grid%ht(grid%itsloc(k),grid%jtsloc(k))
184	#endif
185	END IF
186	#if DM_PARALLEL
187	ts_xlat = wrf_dm_min_real(ts_xlat)
188	ts_xlong = wrf_dm_min_real(ts_xlong)
189	ts_hgt = wrf_dm_min_real(ts_hgt)
190	#endif
191
192	IF ( wrf_dm_on_monitor() ) THEN
193
194	iunit = get_unused_unit()
195	IF ( iunit <= 0 ) THEN
196	CALL wrf_error_fatal('Error in calc_ts_locations: could not find a free Fortran unit.')
197	END IF
198
199	WRITE(grid%ts_filename(k),'(A)') TRIM(grid%nametsloc(grid%id_tsloc(k)))//'.d00.TS'
200	i = LEN_TRIM(grid%ts_filename(k))
201	WRITE(grid%ts_filename(k)(i-4:i-3),'(I2.2)') grid%id
202	OPEN(UNIT=iunit, FILE=TRIM(grid%ts_filename(k)), FORM='FORMATTED', STATUS='REPLACE')
203	#if (EM_CORE == 1)
204	WRITE(UNIT=iunit, &
205	FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2,F6.1,A7)') &
206	grid%desctsloc(grid%id_tsloc(k))//' ', grid%id, grid%id_tsloc(k), &
207	' '//grid%nametsloc(grid%id_tsloc(k)), &
208	' (', grid%lattsloc(grid%id_tsloc(k)), ',', grid%lontsloc(grid%id_tsloc(k)), ') (', &
209	grid%itsloc(k), ',', grid%jtsloc(k), ') (', &
210	ts_xlat, ',', ts_xlong, ') ', &
211	ts_hgt,' meters'
212	#else
213	WRITE(UNIT=iunit, &
214	FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2)') &
215	grid%desctsloc(grid%id_tsloc(k))//' ', grid%id, grid%id_tsloc(k), &
216	' '//grid%nametsloc(grid%id_tsloc(k)), &
217	' (', grid%lattsloc(grid%id_tsloc(k)), ',', grid%lontsloc(grid%id_tsloc(k)), ') (', &
218	grid%itsloc(k), ',', grid%jtsloc(k), ') (', &
219	ts_xlat, ',', ts_xlong, ') '
220	#endif
221	CLOSE(UNIT=iunit)
222	END IF
223	END DO
224
225	END IF
226	#if ((EM_CORE == 1) && (DA_CORE != 1))
227	END IF
228	#endif
229
230	END SUBROUTINE calc_ts_locations
231
232
233	SUBROUTINE calc_ts( grid )
234
235	USE module_domain
236	USE module_model_constants
237
238	IMPLICIT NONE
239
240	! Arguments
241	TYPE (domain), INTENT(INOUT) :: grid
242
243	LOGICAL, EXTERNAL :: wrf_dm_on_monitor
244
245	! Local variables
246	INTEGER :: i, k, mm, n, ix, iy, rc
247	REAL :: earth_u, earth_v, output_t, output_q, clw, xtime_minutes
248	REAL, ALLOCATABLE, DIMENSION(:) :: p8w
249
250	! Parameter ts_model_level:
251	! TRUE to output T, Q, and wind at lowest model level
252	! FALSE to output T and Q at 2-m and wind at 10-m diagnostic levels:
253	LOGICAL, PARAMETER :: ts_model_level = .FALSE.
254
255	IF ( grid%ntsloc_domain .LE. 0 ) RETURN
256
257	#if ((EM_CORE == 1) && (DA_CORE != 1))
258	IF ( grid%dfi_opt /= DFI_NODFI .AND. grid%dfi_stage /= DFI_FST ) RETURN
259	#endif
260
261	n = grid%next_ts_time
262
263	ALLOCATE(p8w(grid%sm32:grid%em32))
264
265	DO i=1,grid%ntsloc_domain
266
267	ix = grid%itsloc(i)
268	iy = grid%jtsloc(i)
269
270	IF (grid%sp31 <= ix .AND. ix <= grid%ep31 .AND. &
271	grid%sp33 <= iy .AND. iy <= grid%ep33) THEN
272
273	IF (ts_model_level) THEN
274
275	!
276	! Output from the lowest model computational level:
277	!
278	#if (EM_CORE == 1)
279	earth_u = grid%u_2(ix,1,iy)grid%cosa(ix,iy)-grid%v_2(ix,1,iy)grid%sina(ix,iy)
280	earth_v = grid%v_2(ix,1,iy)grid%cosa(ix,iy)+grid%u_2(ix,1,iy)grid%sina(ix,iy)
281	output_t = grid%t_2(ix,1,iy)
282	#else
283	earth_u = grid%u(ix,1,iy)
284	earth_v = grid%v(ix,1,iy)
285	output_t = grid%t(ix,1,iy)
286	#endif
287	output_q = grid%moist(ix,1,iy,P_QV)
288
289	ELSE
290
291	!
292	! Output at 2-m and 10-m diagnostic levels:
293	!
294	#if (EM_CORE == 1)
295	earth_u = grid%u10(ix,iy)grid%cosa(ix,iy)-grid%v10(ix,iy)grid%sina(ix,iy)
296	earth_v = grid%v10(ix,iy)grid%cosa(ix,iy)+grid%u10(ix,iy)grid%sina(ix,iy)
297	output_q = grid%q2(ix,iy)
298	#else
299	earth_u = grid%u10(ix,iy)
300	earth_v = grid%v10(ix,iy)
301	output_q = grid%qsfc(ix,iy)
302	#endif
303	output_t = grid%t2(ix,iy)
304
305	END IF
306
307	#if (EM_CORE == 1)
308	! Calculate column-integrated liquid/ice (kg/m^2 or mm)
309	CALL calc_p8w(grid, ix, iy, p8w, grid%sm32, grid%em32)
310	clw=0.
311	DO mm = 1, num_moist
312	IF ( (mm == P_QC) .OR. (mm == P_QR) .OR. (mm == P_QI) .OR. &
313	(mm == P_QS) .OR. (mm == P_QG) ) THEN
314	DO k=grid%sm32,grid%em32-1
315	clw=clw+grid%moist(ix,k,iy,mm)*(p8w(k)-p8w(k+1))
316	END DO
317	END IF
318	END DO
319	clw = clw / g
320	#endif
321
322	CALL domain_clock_get( grid, minutesSinceSimulationStart=xtime_minutes )
323	grid%ts_hour(n,i) = xtime_minutes / 60.
324	grid%ts_u(n,i) = earth_u
325	grid%ts_v(n,i) = earth_v
326	grid%ts_t(n,i) = output_t
327	grid%ts_q(n,i) = output_q
328	grid%ts_psfc(n,i) = grid%psfc(ix,iy)
329	#if (EM_CORE == 1)
330	grid%ts_glw(n,i) = grid%glw(ix,iy)
331	grid%ts_gsw(n,i) = grid%gsw(ix,iy)
332	grid%ts_hfx(n,i) = grid%hfx(ix,iy)
333	grid%ts_lh(n,i) = grid%lh(ix,iy)
334	grid%ts_clw(n,i) = clw
335	grid%ts_rainc(n,i) = grid%rainc(ix,iy)
336	grid%ts_rainnc(n,i) = grid%rainnc(ix,iy)
337	grid%ts_tsk(n,i) = grid%tsk(ix,iy)
338	#else
339	grid%ts_tsk(n,i) = grid%nmm_tsk(ix,iy)
340	#endif
341	grid%ts_tslb(n,i) = grid%tslb(ix,1,iy)
342
343	ELSE
344
345	grid%ts_hour(n,i) = 1.E30
346	grid%ts_u(n,i) = 1.E30
347	grid%ts_v(n,i) = 1.E30
348	grid%ts_t(n,i) = 1.E30
349	grid%ts_q(n,i) = 1.E30
350	grid%ts_psfc(n,i) = 1.E30
351	#if (EM_CORE == 1)
352	grid%ts_glw(n,i) = 1.E30
353	grid%ts_gsw(n,i) = 1.E30
354	grid%ts_hfx(n,i) = 1.E30
355	grid%ts_lh(n,i) = 1.E30
356	grid%ts_clw(n,i) = 1.E30
357	grid%ts_rainc(n,i) = 1.E30
358	grid%ts_rainnc(n,i) = 1.E30
359	#endif
360	grid%ts_tsk(n,i) = 1.E30
361	grid%ts_tslb(n,i) = 1.E30
362
363	END IF
364	END DO
365
366	DEALLOCATE(p8w)
367
368	grid%next_ts_time = grid%next_ts_time + 1
369
370	IF ( grid%next_ts_time > grid%ts_buf_size ) CALL write_ts(grid)
371
372	END SUBROUTINE calc_ts
373
374
375	SUBROUTINE write_ts( grid )
376
377	USE module_domain, ONLY : domain
378	USE module_dm, ONLY : wrf_dm_min_reals
379	USE module_state_description
380
381	IMPLICIT NONE
382
383	! Arguments
384	TYPE (domain), INTENT(INOUT) :: grid
385
386	LOGICAL, EXTERNAL :: wrf_dm_on_monitor
387	INTEGER, EXTERNAL :: get_unused_unit
388
389	! Local variables
390	INTEGER :: i, n, ix, iy, iunit
391	REAL, ALLOCATABLE, DIMENSION(:,:) :: ts_buf
392
393	IF ( grid%ntsloc_domain .LE. 0 ) RETURN
394
395	#if ((EM_CORE == 1) && (DA_CORE != 1))
396	IF ( grid%dfi_opt /= DFI_NODFI .AND. grid%dfi_stage /= DFI_FST ) RETURN
397	#endif
398
399	#ifdef DM_PARALLEL
400	ALLOCATE(ts_buf(grid%ts_buf_size,grid%max_ts_locs))
401
402	ts_buf(:,:) = grid%ts_hour(:,:)
403	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_hour(:,:),grid%ts_buf_size*grid%max_ts_locs)
404
405	ts_buf(:,:) = grid%ts_u(:,:)
406	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_u(:,:),grid%ts_buf_size*grid%max_ts_locs)
407
408	ts_buf(:,:) = grid%ts_v(:,:)
409	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_v(:,:),grid%ts_buf_size*grid%max_ts_locs)
410
411	ts_buf(:,:) = grid%ts_t(:,:)
412	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_t(:,:),grid%ts_buf_size*grid%max_ts_locs)
413
414	ts_buf(:,:) = grid%ts_q(:,:)
415	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_q(:,:),grid%ts_buf_size*grid%max_ts_locs)
416
417	ts_buf(:,:) = grid%ts_psfc(:,:)
418	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_psfc(:,:),grid%ts_buf_size*grid%max_ts_locs)
419
420	#if (EM_CORE == 1)
421	ts_buf(:,:) = grid%ts_glw(:,:)
422	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_glw(:,:),grid%ts_buf_size*grid%max_ts_locs)
423
424	ts_buf(:,:) = grid%ts_gsw(:,:)
425	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_gsw(:,:),grid%ts_buf_size*grid%max_ts_locs)
426
427	ts_buf(:,:) = grid%ts_hfx(:,:)
428	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_hfx(:,:),grid%ts_buf_size*grid%max_ts_locs)
429
430	ts_buf(:,:) = grid%ts_lh(:,:)
431	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_lh(:,:),grid%ts_buf_size*grid%max_ts_locs)
432
433	ts_buf(:,:) = grid%ts_clw(:,:)
434	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_clw(:,:),grid%ts_buf_size*grid%max_ts_locs)
435
436	ts_buf(:,:) = grid%ts_rainc(:,:)
437	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_rainc(:,:),grid%ts_buf_size*grid%max_ts_locs)
438
439	ts_buf(:,:) = grid%ts_rainnc(:,:)
440	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_rainnc(:,:),grid%ts_buf_size*grid%max_ts_locs)
441	#endif
442
443	ts_buf(:,:) = grid%ts_tsk(:,:)
444	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tsk(:,:),grid%ts_buf_size*grid%max_ts_locs)
445
446	ts_buf(:,:) = grid%ts_tslb(:,:)
447	CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tslb(:,:),grid%ts_buf_size*grid%max_ts_locs)
448
449	DEALLOCATE(ts_buf)
450	#endif
451
452	IF ( wrf_dm_on_monitor() ) THEN
453
454	iunit = get_unused_unit()
455	IF ( iunit <= 0 ) THEN
456	CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
457	END IF
458
459	DO i=1,grid%ntsloc_domain
460
461	ix = grid%itsloc(i)
462	iy = grid%jtsloc(i)
463
464	OPEN(UNIT=iunit, FILE=TRIM(grid%ts_filename(i)), STATUS='unknown', POSITION='append', FORM='formatted')
465
466	DO n=1,grid%next_ts_time - 1
467
468	#if (EM_CORE == 1)
469	WRITE(UNIT=iunit,FMT='(i2,f13.6,i5,i5,i5,1x,14(f13.5,1x))') &
470	grid%id, grid%ts_hour(n,i), &
471	grid%id_tsloc(i), ix, iy, &
472	grid%ts_t(n,i), &
473	grid%ts_q(n,i), &
474	grid%ts_u(n,i), &
475	grid%ts_v(n,i), &
476	grid%ts_psfc(n,i), &
477	grid%ts_glw(n,i), &
478	grid%ts_gsw(n,i), &
479	grid%ts_hfx(n,i), &
480	grid%ts_lh(n,i), &
481	grid%ts_tsk(n,i), &
482	grid%ts_tslb(n,i), &
483	grid%ts_rainc(n,i), &
484	grid%ts_rainnc(n,i), &
485	grid%ts_clw(n,i)
486	#else
487	WRITE(UNIT=iunit,FMT='(i2,f13.6,i5,i5,i5,1x,7(f13.5,1x))') &
488	grid%id, grid%ts_hour(n,i), &
489	grid%id_tsloc(i), ix, iy, &
490	grid%ts_t(n,i), &
491	grid%ts_q(n,i), &
492	grid%ts_u(n,i), &
493	grid%ts_v(n,i), &
494	grid%ts_psfc(n,i), &
495	grid%ts_tsk(n,i), &
496	grid%ts_tslb(n,i)
497	#endif
498	END DO
499
500	CLOSE(UNIT=iunit)
501
502	END DO
503
504	END IF
505
506	grid%next_ts_time = 1
507
508	END SUBROUTINE write_ts
509
510
511	#if (EM_CORE == 1)
512	SUBROUTINE calc_p8w(grid, ix, iy, p8w, k_start, k_end)
513
514	USE module_domain
515	USE module_model_constants
516
517	IMPLICIT NONE
518
519	! Arguments
520	TYPE (domain), INTENT(IN) :: grid
521	INTEGER, INTENT(IN) :: ix, iy, k_start, k_end
522	REAL, DIMENSION(k_start:k_end), INTENT(OUT) :: p8w
523
524	! Local variables
525	INTEGER :: k
526	REAL :: z0, z1, z2, w1, w2
527	REAL, DIMENSION(k_start:k_end) :: z_at_w
528	REAL, DIMENSION(k_start:k_end-1) :: z
529
530
531	DO k = k_start, k_end
532	z_at_w(k) = (grid%phb(ix,k,iy)+grid%ph_2(ix,k,iy))/g
533	END DO
534
535	DO k = k_start, k_end-1
536	z(k) = 0.5*(z_at_w(k) + z_at_w(k+1))
537	END DO
538
539	DO k = k_start+1, k_end-1
540	p8w(k) = grid%fnm(k)*(grid%p(ix,k,iy)+grid%pb(ix,k,iy)) + &
541	grid%fnp(k)*(grid%p(ix,k-1,iy)+grid%pb(ix,k-1,iy))
542	END DO
543
544	z0 = z_at_w(k_start)
545	z1 = z(k_start)
546	z2 = z(k_start+1)
547	w1 = (z0 - z2)/(z1 - z2)
548	w2 = 1. - w1
549	p8w(k_start) = w1*(grid%p(ix,k_start,iy)+grid%pb(ix,k_start,iy)) + &
550	w2*(grid%p(ix,k_start+1,iy)+grid%pb(ix,k_start+1,iy))
551
552	z0 = z_at_w(k_end)
553	z1 = z(k_end-1)
554	z2 = z(k_end-2)
555	w1 = (z0 - z2)/(z1 - z2)
556	w2 = 1. - w1
557	p8w(k_end) = exp(w1*log(grid%p(ix,k_end-1,iy)+grid%pb(ix,k_end-1,iy)) + &
558	w2*log(grid%p(ix,k_end-2,iy)+grid%pb(ix,k_end-2,iy)))
559
560	END SUBROUTINE calc_p8w
561	#endif

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format