Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

wrf_timeseries.F @ 3567

Last change on this file since 3567 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

Rev	Line
[2759]	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