Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

wrf_timeseries.F @ 3553

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

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

Download in other formats:

Original Format