Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

iowrf.f @ 23

Last change on this file since 23 was 19, checked in by aslmd, 15 years ago
spiga:mineur
File size: 29.4 KB

Rev	Line
[19]	1	program iowrf
	2	! Program to read/write wrf output.
	3	! OPTION -thin : It will thin data to the ratio given
	4	! OPTION -thina : It will average the fields over a user-specified grid area.
	5	! OPTION -A : De-stagger data
	6	! OPTION -box : Will get data from a user defined box
	7	! Updated: January 8, 2008
	8	! Add de-staggering option
	9	! Updated: Jun 19, 2007
	10	! Change time to unlimted on output
	11	! Origincal code:
	12	! Cindy Bruyere - March 2006
	13	! Some code borrowed from Jim Bresch
	14	!=================================Run Program================================
	15	! To extract a box from your input file
	16	! iowrf wrfout_file -box x 10 50 y 10 60 [-debug]
	17	!
	18	! To thin your input file down (by picking up corresponding points)
	19	! iowrf wrfout_file -thin 3 [-debug]
	20	!
	21	! To thin your input file down (by averaging)
	22	! iowrf wrfout_file -thina 3 [-debug]
	23	!
	24	! To de-stagger the data
	25	! iowrf wrfout_file -A
	26	!
	27	! To CREATE large 64bit data files
	28	! -64bit
	29	!
	30	! To see more options
	31	! iowrf -help
	32	!
	33	! The output will be written to a file with original file name + -box/-thin/-thina/-A
	34	!
	35	!=================================Make Executable============================
	36	! Make executable:
	37	! DEC Alpha
	38	! f90 iowrf.f -L/usr/local/netcdf/lib -lnetcdf -lm \
	39	! -I/usr/local/netcdf/include -free -o iowrf
	40	!
	41	! linux flags
	42	! pgf90 iowrf.f -L/usr/local/netcdf/lib -lnetcdf -lm \
	43	! -I/usr/local/netcdf/include -Mfree -o iowrf
	44	!
	45	! Sun flags
	46	! f90 iowrf.f -L/usr/local/netcdf/lib -lnetcdf -lm \
	47	! -I/usr/local/netcdf/include -free -o iowrf
	48	!
	49	! SGI flags
	50	! f90 iowrf.f -L/usr/local/netcdf/lib -lnetcdf -lm \
	51	! -I/usr/local/netcdf/include -freeform -o iowrf
	52	!
	53	! IBM flags (NCAR bluesky - 32bit)
	54	! xlf iowrf.f -L/usr/local/lib32/r4i4 -lnetcdf -lm \
	55	! -I/usr/local/include -qfree=f90 -o iowrf
	56	!
	57	! IBM flags (NCAR bluesky - 64bit)
	58	! xlf iowrf.f -L/usr/local/lib64/r4i4 -lnetcdf -lm \
	59	! -I/usr/local/include -qfree=f90 -o iowrf
	60	!
	61	! IBM flags (NCAR bluevista)
	62	! xlf iowrf.f -L/usr/local/netcdf/lib -lnetcdf -lm \
	63	! -I/usr/local/netcdf/include -qfree=f90 -o iowrf
	64	!
	65	! Mac flags (with xlf compiler)
	66	! xlf iowrf.f -L/usr/local/netcdf-xlf/lib -lnetcdf -lm \
	67	! -I/usr/local/netcdf-xlf/include -qfree=f90 -o iowrf
	68	!
	69	! Mac flags (with g95 compiler)
	70	! g95 iowrf.f -L/usr/local/netcdf/lib -lnetcdf -lm \
	71	! -I/usr/local/netcdf/include -ffree-form -o iowrf
	72	!
	73	!============================================================================
	74
	75	implicit none
	76	INCLUDE 'netcdf.inc'
	77	integer :: jdim
	78	parameter (jdim=6)
	79	integer ncid, status
	80	integer ishape(jdim)
	81	integer ishape2(jdim)
	82	character cval*50
	83	character name*31
	84	character (len=31),allocatable, dimension(:) :: dname
	85	integer, allocatable, dimension(:) :: dval, dval2
	86	real, allocatable, dimension(:,:,:,:) :: data, data2
	87	double precision, allocatable, dimension(:,:,:,:) :: ddata, ddata2
	88	integer, allocatable, dimension(:,:,:,:) :: idata, idata2
	89	character, allocatable, dimension(:,:,:,:) :: text
	90	character omit(10)*80
	91	integer :: start_dims(4)
	92	integer :: dims_in(4), dims_out(4), box_start(3), box_end(3)
	93	integer :: firstS,firstE, secondS,secondE, thirdS,thirdE
	94	integer :: idm, ndims, nvars, natt, ngatts, nunlimdimid, iratio
	95	integer :: i, ii, j, iweg, jweg, isng, jsng, ibtg, jbtg, ix, iy
	96	integer :: i_shape_we, i_shape_sn, i_shape_bt
	97	integer :: i_shape_we_stag, i_shape_sn_stag, i_shape_bt_stag
	98	integer :: ilen, itype, ival, na
	99	integer :: mcid
	100	real :: dx, rval
	101	real :: new_cen
	102	real :: okm
	103	character (len=80) :: input_file, output_file
	104	character (len=10) :: option
	105	logical :: debug=.FALSE.
	106	logical :: x_ave=.FALSE.
	107	logical :: y_ave=.FALSE.
	108	logical :: bit64
	109
	110	call read_args(input_file,option,iratio,box_start,box_end,bit64,debug)
	111	output_file = trim(input_file)//option
	112
	113	write(6,*)
	114	write(6,*) "#########################################"
	115	write(6,*) "Running IOWRF "
	116	write(6,*)
	117	write(6,*) "INPUT FILE: ",trim(input_file)
	118	write(6,*) "OUTPUT FILE: ",trim(output_file)
	119	write(6,*) "OPTION: ",option
	120
	121	IF (debug) THEN
	122	if ( option(1:5) == '-thin' ) then ! used for -thina and -thin
	123	write(6,*) "RATIO: ",iratio
	124	elseif ( option == '-box' )then
	125	write(6,*) "BOX START (x y z): ",box_start
	126	write(6,*) "BOX END (x y z): ",box_end
	127	endif
	128	ENDIF
	129
	130
	131	! OPEN INPUT AND OUTPUT FILE
	132	! output_file is input_file_new
	133	status = nf_open(input_file, 0, ncid)
	134	if (status .ne. nf_noerr) call handle_err(status)
	135	if (bit64) then
	136	status = nf_create(output_file, NF_64BIT_OFFSET, mcid)
	137	else
	138	status = nf_create(output_file, 0, mcid)
	139	endif
	140	if (status .ne. nf_noerr) call handle_err(status)
	141
	142	! GET BASIC INFORMTION ABOUT THE FILE
	143	! most important
	144	! ndims: number of dimensions
	145	! nvars: number of variables
	146	! ngatts: number of global attributes
	147	status = nf_inq(ncid, ndims, nvars, ngatts, nunlimdimid)
	148	if (status .ne. nf_noerr) call handle_err(status)
	149	IF (debug) THEN
	150	write(6,*)
	151	write(6,'(4(A,i4))') ' ndims = ',ndims,' nvars = ',nvars,' ngatts = ',ngatts, &
	152	' nunlimdimid =',nunlimdimid
	153	write(6,*)
	154	ENDIF
	155
	156	! ALLOCATE SOME VARIABLES
	157	allocate (dval(ndims))
	158	allocate(dval2(ndims))
	159	allocate(dname(ndims))
	160
	161	! GET SOME BASIC DIMS FROM INPUT_FILE
	162	dx = -99.
	163	status = nf_get_att_real (ncid, nf_global, 'DX', dx)
	164	status = nf_get_att_int (ncid, nf_global, 'WEST-EAST_GRID_DIMENSION', iweg)
	165	status = nf_get_att_int (ncid, nf_global, 'SOUTH-NORTH_GRID_DIMENSION', isng)
	166	status = nf_get_att_int (ncid, nf_global, 'BOTTOM-TOP_GRID_DIMENSION', ibtg)
	167	IF (debug) THEN
	168	write(6,*) "BASICS from input file:"
	169	write(6,*) " DX= ", dx
	170	write(6,*) " X= ", iweg
	171	write(6,*) " Y= ", isng
	172	write(6,*) " Z= ", ibtg
	173	ENDIF
	174	if (dx .lt. 0.) stop 'dx is bad'
	175
	176	! CALCULATE DIMS FOR OUTPUT FILE
	177	IF ( option(1:5) == '-thin' ) THEN ! used for -thina and -thin
	178	okm = dx*iratio
	179	jweg = int((iweg-1)/iratio) + 1
	180	jsng = int((isng-1)/iratio) + 1
	181	jbtg = ibtg
	182	ELSEIF ( option == '-box' ) THEN
	183	okm = dx
	184	jweg = iweg
	185	jsng = isng
	186	jbtg = ibtg
	187	if ( box_end(1) .ne. 0 ) jweg = int(box_end(1) - box_start(1)) + 1
	188	if ( box_end(2) .ne. 0 ) jsng = int(box_end(2) - box_start(2)) + 1
	189	if ( box_end(3) .ne. 0 ) jbtg = int(box_end(3) - box_start(3)) + 1
	190	ELSE
	191	okm = dx
	192	jweg = iweg
	193	jsng = isng
	194	jbtg = ibtg
	195	ENDIF
	196	IF (debug) THEN
	197	write(6,*) "BASICS for output file:"
	198	write(6,*) " DX= ", okm
	199	write(6,*) " X= ", jweg
	200	write(6,*) " Y= ", jsng
	201	write(6,*) " Z= ", jbtg
	202	ENDIF
	203	!! We also need to fix the CEN_LAT and CEN_LON later, so get
	204	!! the middle of the new domain
	205	ix = int((jweg-1)/2.)
	206	iy = int((jsng-1)/2.)
	207	if ( ix .eq. int(jweg/2.) ) x_ave = .TRUE.
	208	if ( iy .eq. int(jsng/2.) ) y_ave = .TRUE.
	209	ix = int(jweg/2.)
	210	iy = int(jsng/2.)
	211
	212	! READ ALL DIMS FROM INPUT FILE AND CREATE DIMS FOR OUTPUT FILE
	213	IF (debug) THEN
	214	write(6,*)
	215	write(6,*) "FILE dimensions:"
	216	ENDIF
	217	i_shape_we = 0
	218	i_shape_sn = 0
	219	i_shape_bt = 0
	220	i_shape_we_stag = 0
	221	i_shape_sn_stag = 0
	222	i_shape_bt_stag = 0
	223
	224	do i = 1, ndims
	225	status = nf_inq_dim(ncid, i, dname(i), dval(i))
	226	dval2(i) = dval(i)
	227	! CAUTION -- this stuff is hard-wired
	228	if (dname(i) .eq. 'west_east_stag') then
	229	dval2(i) = jweg
	230	i_shape_we_stag = i
	231	else if (dname(i) .eq. 'west_east') then
	232	dval2(i) = jweg-1
	233	i_shape_we = i
	234	else if (dname(i) .eq. 'south_north_stag') then
	235	dval2(i) = jsng
	236	i_shape_sn_stag = i
	237	else if (dname(i) .eq. 'south_north') then
	238	dval2(i) = jsng-1
	239	i_shape_sn = i
	240	else if (dname(i) .eq. 'bottom_top_stag') then
	241	dval2(i) = jbtg
	242	i_shape_bt_stag = i
	243	else if (dname(i) .eq. 'bottom_top') then
	244	dval2(i) = jbtg-1
	245	i_shape_bt = i
	246	endif
	247	if ( dname(i) == "Time" ) then
	248	status = nf_def_dim(mcid, dname(i), NF_UNLIMITED, i)
	249	else
	250	status = nf_def_dim(mcid, dname(i), dval2(i), i)
	251	end if
	252	IF (debug) THEN
	253	write(6,'(i4," : ",A," in = ",i4," (out = ",i4,")")') &
	254	i,dname(i),dval(i),dval2(i)
	255	ENDIF
	256	enddo
	257	IF (.not. debug) THEN
	258	write(6,*)
	259	write(6,*) "Set up file DIMENSIONS"
	260	ENDIF
	261
	262	! DEALING WITH THE GLOBAL ATTRIBUTES
	263	IF (debug) THEN
	264	write(6,*)
	265	write(6,*) "FILE attributes:"
	266	ENDIF
	267	do i = 1, ngatts
	268	status = nf_inq_attname(ncid, nf_global, i, name)
	269	status = nf_inq_atttype(ncid, nf_global, name, itype)
	270	status = nf_inq_attlen(ncid, nf_global, name, ilen)
	271
	272	if ( itype .eq. 2 ) then ! characters
	273	status = nf_get_att_text (ncid, nf_global, name, cval)
	274	IF (debug) THEN
	275	write(6,'(i4," : ",A," in = ",A," (out = ",$)') &
	276	i,name,cval(1:ilen)
	277	ENDIF
	278	if(name(1:5) .eq. 'TITLE') then
	279	cval = cval(1:ilen)//" : iowrf"//option
	280	ilen = len_trim(cval)
	281	endif
	282	IF (debug) write(6,'(A,")")') cval(1:ilen)
	283	status = nf_put_att_text(mcid, nf_global, name, ilen,&
	284	cval(1:ilen))
	285
	286	elseif ( itype .eq. 4 ) then ! integers
	287	status = nf_get_att_int (ncid, nf_global, name, ival)
	288	IF (debug) THEN
	289	write(6,'(i4," : ",A," in = ",i7," (out = ",$)') &
	290	i,name,ival
	291	ENDIF
	292	if(name .eq. 'WEST-EAST_GRID_DIMENSION') ival = jweg
	293	if(name .eq. 'SOUTH-NORTH_GRID_DIMENSION') ival = jsng
	294	if(name .eq. 'BOTTOM-TOP_GRID_DIMENSION') ival = jbtg
	295	IF (debug) write(6,'(i7,")")') ival
	296	status = nf_put_att_int(mcid, nf_global, name, itype,&
	297	ilen, ival)
	298
	299	elseif ( itype .eq. 5 ) then ! real
	300	status = nf_get_att_real (ncid, nf_global, name, rval)
	301	IF (debug) THEN
	302	write(6,'(i4," : ",A," in = ",G18.10E2," (out = ",$)') &
	303	i,name,rval
	304	ENDIF
	305	if(name(1:2) .eq. 'DX' .or. name(1:2) .eq. 'DY') rval = okm
	306	IF (debug) write(6,'(G18.10E2,")")') rval
	307	status = nf_put_att_real(mcid, nf_global, name, itype,&
	308	ilen, rval)
	309	endif
	310	enddo
	311	IF ( .not. debug ) THEN
	312	write(6,*) "Write file ATTRIBUTES"
	313	write(6,*)
	314	ENDIF
	315
	316
	317	! TRAIN FILE
	318	do i = 1, nvars
	319	status = nf_inq_var(ncid, i, cval, itype, idm, ishape, natt)
	320	ishape2 = ishape
	321	if ( idm .ge. 4 ) then
	322	do ii=1,idm
	323	IF ( option == "-A" .AND. ishape2(ii) == i_shape_bt_stag ) THEN
	324	ishape2(ii) = i_shape_bt
	325	ELSEIF ( option == "-A" .AND. ishape2(ii) == i_shape_we_stag ) THEN
	326	ishape2(ii) = i_shape_we
	327	ELSEIF ( option == "-A" .AND. ishape2(ii) == i_shape_sn_stag ) THEN
	328	ishape2(ii) = i_shape_sn
	329	END IF
	330	enddo
	331	end if
	332
	333	status = nf_def_var(mcid, cval, itype, idm, ishape2, i)
	334	do na = 1, natt
	335	status = nf_inq_attname(ncid, i, na, name)
	336	status = nf_copy_att(ncid, i, name, mcid, i)
	337	enddo
	338	enddo
	339	status = nf_enddef(mcid)
	340
	341	! ########## LOOP THROUGH THE DATA
	342	IF (debug) THEN
	343	write(6,*)
	344	write(6,*)
	345	ENDIF
	346	write(6,*) "Write file VARIABLES:"
	347	start_dims = 1
	348	do i = 1, nvars
	349	status = nf_inq_var(ncid, i, cval, itype, idm, ishape, natt)
	350	ishape2 = ishape
	351	if ( idm .ge. 4 ) then
	352	do ii=1,idm
	353	IF ( option == "-A" .AND. ishape2(ii) == i_shape_bt_stag ) THEN
	354	ishape2(ii) = i_shape_bt
	355	ELSEIF ( option == "-A" .AND. ishape2(ii) == i_shape_we_stag ) THEN
	356	ishape2(ii) = i_shape_we
	357	ELSEIF ( option == "-A" .AND. ishape2(ii) == i_shape_sn_stag ) THEN
	358	ishape2(ii) = i_shape_sn
	359	END IF
	360	enddo
	361	end if
	362	IF (debug) THEN
	363	write(6,*)
	364	ENDIF
	365	write(6,*) 'VARIABLE: ',trim(cval)
	366
	367	! GET THE DIMS FOR INPUT AND OUTPUT FROM THE SHAPE
	368	dims_in = 1
	369	dims_out = 1
	370	do ii = 1,idm
	371	dims_in(ii) = dval(ishape(ii))
	372	dims_out(ii) = dval2(ishape2(ii))
	373	enddo
	374	IF (debug) THEN
	375	write(6,*) ' DIMS IN: ',dims_in
	376	write(6,*) ' DIMS OUT: ',dims_out
	377	ENDIF
	378
	379	IF ( option == '-box' ) THEN
	380	!! Get the start and end dimensions of the box in the input file
	381	firstS = 1
	382	firstE = dims_out(1)
	383	secondS = 1
	384	secondE = dims_out(2)
	385	thirdS = 1
	386	thirdE = dims_out(3)
	387
	388	if (idm.eq.2 .and. dims_out(1).ge.jbtg-1 .and. box_end(3).ne.0) then
	389	firstS = box_start(3)
	390	firstE = box_end(3)
	391	if (dims_out(3) .eq. jbtg-1) firstE = firstE-1
	392	endif
	393	if (idm .ge. 3) then
	394	if (box_end(1) .ne. 0) then
	395	firstS = box_start(1)
	396	firstE = box_end(1)
	397	if (dims_out(1) .eq. jweg-1) firstE = firstE-1
	398	endif
	399	if (box_end(2) .ne. 0) then
	400	secondS = box_start(2)
	401	secondE = box_end(2)
	402	if (dims_out(2) .eq. jsng-1) secondE = secondE-1
	403	endif
	404	if (idm == 4 .and. box_end(3).ne.0) then
	405	thirdS = box_start(3)
	406	thirdE = box_end(3)
	407	if (dims_out(3) .eq. jbtg-1) thirdE = thirdE-1
	408	endif
	409	endif
	410	ENDIF
	411
	412	! ALLOCATE THE INPUT AND OUTPUT ARRAYS
	413	! READ THE DATA FROM INPUT FILE
	414	! THIN THE GRID IF NEEDED, OR GET THE CORRECT BOX
	415
	416	IF (itype .eq. 2) THEN ! character
	417	allocate (text(dims_in(1), dims_in(2), dims_in(3), &
	418	dims_in(4)))
	419	status = nf_get_var_text(ncid, i, text)
	420	IF (debug) write(6,*) ' SAMPLE VALUE = ',text(:,1,1,1)
	421	status = nf_put_vara_text (mcid, i, start_dims, dims_in, text)
	422	deallocate (text)
	423
	424	ELSEIF (itype .eq. 4) THEN ! integer
	425	allocate (idata(dims_in(1), dims_in(2), dims_in(3), &
	426	dims_in(4)))
	427	allocate(idata2(dims_out(1),dims_out(2),dims_out(3),&
	428	dims_out(4)))
	429	status = nf_get_var_int(ncid, i, idata)
	430	IF (debug) write(6,*) ' SAMPLE VALUE = ',idata(1,1,1,1)
	431	IF ( option == '-thina' ) THEN
	432	if (idm .ge. 3) then
	433	IF (debug) write(6,*) ' Grid is thinned with a ratio of ',iratio
	434	allocate(data2(dims_out(1),dims_out(2),dims_out(3),&
	435	dims_out(4)))
	436	call thin_ave (real(idata),data2,dims_in(1),dims_in(2), &
	437	dims_in(3),dims_in(4),dims_out(1),dims_out(2), &
	438	dims_out(3),dims_out(4),iratio)
	439	idata2 = int(data2)
	440	deallocate(data2)
	441	else
	442	idata2 = idata
	443	endif
	444	ELSEIF ( option == '-thin' ) THEN
	445	if (idm .ge. 3) then
	446	IF (debug) write(6,*) ' Grid is thinned with a ratio of ',iratio
	447	idata2 = idata(1:dims_in(1):iratio,1:dims_in(2):iratio,:,:)
	448	else
	449	idata2 = idata
	450	endif
	451	ELSEIF ( option == '-A' ) THEN
	452	idata2 = idata
	453	ELSEIF ( option == '-box') THEN
	454	IF (debug) write(6,*) ' a BOX is extracted from the input domain '
	455	idata2 = idata(firstS:firstE,secondS:secondE,thirdS:thirdE,:)
	456	ENDIF
	457	status = nf_put_vara_int (mcid, i, start_dims, dims_out, idata2)
	458	deallocate (idata)
	459	deallocate (idata2)
	460
	461	ELSEIF (itype .eq. 5) THEN ! real
	462	allocate (data(dims_in(1), dims_in(2), dims_in(3), &
	463	dims_in(4)))
	464	allocate(data2(dims_out(1),dims_out(2),dims_out(3), &
	465	dims_out(4)))
	466	status = nf_get_var_real(ncid, i, data)
	467	IF (debug) write(6,*) ' SAMPLE VALUE = ',data(1,1,1,1)
	468	IF ( option == '-thina' ) THEN
	469	if (idm .ge. 3) then
	470	IF (debug) write(6,*) ' Grid is thinned with a ratio of ',iratio
	471	call thin_ave (data,data2,dims_in(1),dims_in(2), &
	472	dims_in(3),dims_in(4),dims_out(1),dims_out(2), &
	473	dims_out(3),dims_out(4),iratio)
	474	else
	475	data2 = data
	476	endif
	477	ELSEIF ( option == '-thin' ) THEN
	478	if (idm .ge. 3) then
	479	IF (debug) write(6,*) ' Grid is thinned with a ratio of ',iratio
	480	data2 = data(1:dims_in(1):iratio,1:dims_in(2):iratio,:,:)
	481	else
	482	data2 = data
	483	endif
	484	ELSEIF ( option == '-A' ) THEN
	485	if (idm .eq. 4 .AND. (dims_in(1) > dims_out(1)) ) then
	486	IF (debug) write(6,*) ' de-staggering in the X direction'
	487	data2 = (data(1:dims_in(1)-1,:,:,:)+data(2:dims_in(1),:,:,:))*.5
	488	elseif (idm .eq. 4 .AND. (dims_in(2) > dims_out(2)) ) then
	489	IF (debug) write(6,*) ' de-staggering in the Y direction'
	490	data2 = (data(:,1:dims_in(2)-1,:,:)+data(:,2:dims_in(2),:,:))*.5
	491	elseif (idm .eq. 4 .AND. (dims_in(3) > dims_out(3)) ) then
	492	IF (debug) write(6,*) ' de-staggering in the Y direction'
	493	data2 = (data(:,:,1:dims_in(3)-1,:)+data(:,:,3:dims_in(2),:))*.5
	494	else
	495	data2 = data
	496	endif
	497	ELSEIF ( option == '-box') THEN
	498	IF (debug) write(6,*) ' a BOX is extracted from the input domain '
	499	data2 = data(firstS:firstE,secondS:secondE,thirdS:thirdE,:)
	500	ENDIF
	501	status = nf_put_vara_real (mcid, i, start_dims, dims_out, data2)
	502	IF ( cval == 'XLAT' .or. cval == 'XLONG' ) THEN
	503	! We need fix the box's center long and lat
	504	new_cen = data2(ix,iy,1,1)
	505	if ( x_ave .and. y_ave ) then
	506	new_cen = (data2(ix, iy,1,1)+data2(ix ,iy+1,1,1)+ &
	507	data2(ix+1,iy,1,1)+data2(ix+1,iy+1,1,1))/4.
	508	elseif ( x_ave .and. .not. y_ave ) then
	509	new_cen = (data2(ix, iy,1,1)+data2(ix+1,iy ,1,1))/2.
	510	elseif ( .not. x_ave .and. y_ave ) then
	511	new_cen = (data2(ix, iy,1,1)+data2(ix ,iy+1,1,1))/2.
	512	endif
	513	ENDIF
	514	IF ( cval == 'XLAT' ) THEN
	515	IF (debug) write(6,*) ' Fix global attribute CEN_LAT: now = ', new_cen
	516	status = nf_inq_atttype(ncid, nf_global, 'CEN_LAT', itype)
	517	status = nf_inq_attlen(ncid, nf_global, 'CEN_LAT', ilen)
	518	status = nf_put_att_real(mcid, nf_global, 'CEN_LAT', itype,&
	519	ilen, new_cen)
	520	ELSEIF ( cval == 'XLONG' ) THEN
	521	IF (debug) write(6,*) ' Fix global attribute CEN_LON: now = ', new_cen
	522	status = nf_inq_atttype(ncid, nf_global, 'CEN_LON', itype)
	523	status = nf_inq_attlen(ncid, nf_global, 'CEN_LON', ilen)
	524	status = nf_put_att_real(mcid, nf_global, 'CEN_LON', itype,&
	525	ilen, new_cen)
	526	ENDIF
	527	deallocate (data)
	528	deallocate (data2)
	529
	530	ELSEIF (itype .eq. 6) THEN ! double
	531	allocate (ddata(dims_in(1), dims_in(2), dims_in(3), &
	532	dims_in(4)))
	533	allocate(ddata2(dims_out(1),dims_out(2),dims_out(3),&
	534	dims_out(4)))
	535	status = nf_get_var_double(ncid, i, ddata)
	536	IF (debug) write(6,*) ' SAMPLE VALUE = ',ddata(1,1,1,1)
	537	IF ( option == '-thina' ) THEN
	538	if (idm .ge. 3) then
	539	IF (debug) write(6,*) ' Grid is thinned with a ratio of ',iratio
	540	allocate(data2(dims_out(1),dims_out(2),dims_out(3),&
	541	dims_out(4)))
	542	call thin_ave (real(ddata),data2,dims_in(1),dims_in(2), &
	543	dims_in(3),dims_in(4),dims_out(1),dims_out(2), &
	544	dims_out(3),dims_out(4),iratio)
	545	ddata2 = data2
	546	deallocate (data2)
	547	else
	548	ddata2 = ddata
	549	endif
	550	ELSEIF ( option == '-thin' ) THEN
	551	if (idm .ge. 3) then
	552	IF (debug) write(6,*) ' Grid is thinned with a ratio of ',iratio
	553	ddata2 = ddata(1:dims_in(1):iratio,1:dims_in(2):iratio,:,:)
	554	else
	555	ddata2 = ddata
	556	endif
	557	ELSEIF ( option == '-A' ) THEN
	558	if (idm .eq. 4 .AND. (dims_in(1) > dims_out(1)) ) then
	559	IF (debug) write(6,*) ' de-staggering in the X direction'
	560	ddata2 = (ddata(1:dims_in(1)-1,:,:,:)+ddata(2:dims_in(1),:,:,:))*.5
	561	elseif (idm .eq. 4 .AND. (dims_in(2) > dims_out(2)) ) then
	562	IF (debug) write(6,*) ' de-staggering in the Y direction'
	563	ddata2 = (ddata(:,1:dims_in(2)-1,:,:)+ddata(:,2:dims_in(2),:,:))*.5
	564	elseif (idm .eq. 4 .AND. (dims_in(3) > dims_out(3)) ) then
	565	IF (debug) write(6,*) ' de-staggering in the Z direction'
	566	ddata2 = (ddata(:,:,1:dims_in(3)-1,:)+ddata(:,:,3:dims_in(2),:))*.5
	567	else
	568	ddata2 = ddata
	569	endif
	570	ELSEIF ( option == '-box') THEN
	571	IF (debug) write(6,*) ' a BOX is extracted from the input domain '
	572	ddata2 = ddata(firstS:firstE,secondS:secondE,thirdS:thirdE,:)
	573	ENDIF
	574	status = nf_put_vara_double (mcid, i, start_dims, dims_out, ddata2)
	575	deallocate (ddata)
	576	deallocate (ddata2)
	577	ELSE
	578	stop 'trouble - do not know the variable type'
	579	ENDIF
	580
	581	ENDDO ! END OF VARIABLE LOOP
	582	status = nf_close(mcid)
	583
	584	write(6,*)
	585	write(6,*) "SUCCESS - we are out of here"
	586	write(6,*) "#########################################"
	587
	588	end program iowrf
	589	!---------------------------------------------------------------------
	590	subroutine handle_err(status)
	591	integer status
	592	write(6,*) 'Error number ',status
	593	stop
	594	end subroutine
	595	!---------------------------------------------------------------------
	596	subroutine thin_ave (ain, aou, a1, a2, a3, a4, b1, b2, b3, b4, &
	597	iratio)
	598	! average one array into another in x,y.
	599	integer a1, a2, a3, a4, b1, b2, b3, b4,iratio
	600	real ain(a1,a2,a3,a4), aou(b1,b2,b3,b4)
	601	!write(6,*) 'begin thin_ave, ratio = ',iratio
	602	!write(6,*) 'a1 = ',a1,' a2 = ',a2,' a3 = ',a3,' a4 = ',a4
	603	!write(6,*) 'b1 = ',b1,' b2 = ',b2,' b3 = ',b3,' b4 = ',b4
	604	do k4 = 1, b4
	605	do k3 = 1, b3
	606	do j = 1, b2
	607	ymx = (((j-1) * iratio) + 1 ) + iratio/2.
	608	ymn = (((j-1) * iratio) + 1 ) - iratio/2.
	609	ymx = amin1(float(a2),ymx)
	610	ymn = amax1(1.,ymn)
	611	!write(6,*) 'ymn = ',ymn,' ymx = ',ymx
	612	do i = 1, b1
	613	!write(6,*) 'i = ',i,' j = ',j
	614	xmx = (((i-1) * iratio) + 1 ) + iratio/2.
	615	xmn = (((i-1) * iratio) + 1 ) - iratio/2.
	616	xmx = amin1(float(a1),xmx)
	617	xmn = amax1(1.,xmn)
	618	!write(6,*) 'xmn = ',xmn,' xmx = ',xmx
	619	nc = 0
	620	sum = 0.
	621	nn1 = int(ymn+.5)
	622	nn2 = int(ymx)
	623	do n = nn1, nn2
	624	mm1 = int(xmn+.5)
	625	mm2 = int(xmx)
	626	!write(6,*) 'nn1 = ',nn1,' nn2 = ',nn2
	627	!write(6,*) 'mm1 = ',mm1,' mm2 = ',mm2
	628	do m = mm1, mm2
	629	sum = ain(m,n,k3,k4) + sum
	630	!write(6,*) 'm = ',m,' n = ',n,' ain = ',ain(m,n,k3,k4),&
	631	! ' sum = ',sum
	632	nc = nc + 1
	633	enddo
	634	enddo
	635	aou(i,j,k3,k4) = sum/float(nc)
	636	!write(6,*) 'i = ',i,' value = ',aou(i,j,k3,k4)
	637	enddo
	638	enddo
	639	enddo
	640	enddo
	641
	642	end subroutine thin_ave
	643	!--------------------------------------------------------
	644
	645	subroutine read_args(input_file,option,iratio,box_start,box_end,bit64,debug)
	646
	647	implicit none
	648	character (len=80) :: input_file
	649	character (len=10) :: option
	650
	651	integer :: iratio
	652	integer :: box_start(3), box_end(3)
	653	logical :: bit64, debug
	654
	655	integer :: numarg, i, idummy, idummy1, idummy2
	656	real :: rdummy
	657	integer, external :: iargc
	658	character (len=80) :: dummy, dir
	659
	660	! set up some defaults first
	661	input_file = " "
	662	option = " "
	663	idummy1 = 0
	664	idummy2 = 0
	665	box_start = 0
	666	box_end = 0
	667	bit64 = .FALSE.
	668	numarg = iargc()
	669	i = 1
	670
	671	if (numarg .lt. 1) call help_info
	672
	673	do while (i <= numarg)
	674	call getarg(i,dummy)
	675
	676	if (dummy(1:1) == "-") then ! We have an option, else it is the filename
	677
	678	SELECTCASE (trim(dummy))
	679	CASE ("-help")
	680	call help_info
	681	CASE ("-h")
	682	call help_info
	683	CASE ("-debug")
	684	debug = .TRUE.
	685	CASE ("-thina")
	686	option = dummy
	687	i = i+1
	688	call getarg(i,dummy)
	689	read(dummy,'(i3)')idummy
	690	iratio = idummy
	691	if ( iratio .lt. 2 ) STOP ' Must supply a ratio of 2 or more '
	692	CASE ("-thin")
	693	option = dummy
	694	i = i+1
	695	call getarg(i,dummy)
	696	read(dummy,'(i3)')idummy
	697	iratio = idummy
	698	if ( iratio .lt. 2 ) STOP ' Must supply a ratio of 2 or more '
	699	CASE ("-A")
	700	option = dummy
	701	CASE ("-box")
	702	option = dummy
	703	DO
	704	i = i+1
	705	call getarg(i,dir)
	706	if (dir(1:1) == '-') then
	707	i=i-1
	708	exit
	709	endif
	710	if (dir.ne.'x') then
	711	if (dir.ne.'y') then
	712	if (dir.ne.'z') then
	713	i=i-1
	714	exit
	715	endif
	716	endif
	717	endif
	718	i = i+1
	719	call getarg(i,dummy)
	720	read(dummy,'(i3)')idummy1
	721	i = i+1
	722	call getarg(i,dummy)
	723	read(dummy,'(i3)')idummy2
	724	if (dir.eq.' '.or.idummy1.eq.0.or.idummy2.eq.0) exit
	725	if ( dir == 'x' ) then
	726	box_start(1) = idummy1
	727	box_end(1) = idummy2
	728	endif
	729	if ( dir == 'y' ) then
	730	box_start(2) = idummy1
	731	box_end(2) = idummy2
	732	endif
	733	if ( dir == 'z' ) then
	734	box_start(3) = idummy1
	735	box_end(3) = idummy2
	736	endif
	737	idummy1 = 0
	738	idummy2 = 0
	739	ENDDO
	740	CASE ("-64bit")
	741	bit64 = .TRUE.
	742	CASE DEFAULT
	743	call help_info
	744	END SELECT
	745	else
	746	input_file = dummy
	747	endif
	748
	749	i = i+1
	750
	751	enddo
	752
	753	if (input_file == " ") call help_info
	754
	755	end subroutine read_args
	756	!------------------------------------------------------------------------------
	757
	758	subroutine help_info
	759
	760	print*," "
	761	print*," iowrf wrf_data_file_name [-options] "
	762	print*," "
	763	print*," Current options available are:"
	764	print*," -help : Print this information"
	765	print*," -h : Print this information"
	766	print*," "
	767	print*," -thina x : Thin the input grid, with a ratio of x"
	768	print*," Example:"
	769	print*," -thina 3 "
	770	print*," Will thin the grid with a ratio of 3, i.e., "
	771	print*," a 12km grid will be upscaled to a 36km grid"
	772	print*," The new grid point will be an average of the surrounding points"
	773	print*," "
	774	print*," -thin x : Thin the input grid, with a ratio of x"
	775	print*," Example:"
	776	print*," -thin 3 "
	777	print*," Will thin the grid with a ratio of 3, i.e., "
	778	print*," a 12km grid will be upscaled to a 36km grid"
	779	print*," The new grid point will be the feedback value from the input domain"
	780	print*," "
	781	print*," -box [ ] : Will extract a box out of the input grid"
	782	print*," The box can have values for x/y/z "
	783	print*," Examples:"
	784	print*," -box x 10 30 y 20 40 z 5 15"
	785	print*," -box x 10 30 y 20 40 "
	786	print*," -box x 10 30 "
	787	print*," -box y 20 40 "
	788	print*," "
	789	print*," -A : De-stagger output"
	790	print*," "
	791	print*," -64bit : To create large 64bit data files"
	792	print*," "
	793	end subroutine help_info
	794
	795	!------------------------------------------------------------------------------
	796

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

Context Navigation

source: trunk/mesoscale/LMD_MM_MARS/SRC/POSTPROC/iowrf.f @ 23

Download in other formats: