Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

hrecast.F90 @ 3436

Last change on this file since 3436 was 1905, checked in by mturbet, 7 years ago
import some generic utilities from Mars GCM
File size: 48.0 KB

Rev	Line
[1905]	1
	2
	3	program hrecast
	4
	5	! This program reads fields from GCM output files
	6	! (diagfi.nc, stats.nc, concat.nc, ...)
	7	! and recasts it on a new horizontal grid specified by the user.
	8	! The output file name is automatically generated from input file name:
	9	! if input file name is 'input.nc', then output file will be 'input_h.nc'
	10	!
	11	! NB: It is OK if the output grid is not like typical GCM outputs grids,
	12	! but longitudes range must be something like from -180 to 180 (although
	13	! you may omit either of these endpoints), and latitudes must range
	14	! from 90 to -90 (again, you may omit the endpoints).
	15	! EM 09/2009
	16	! TN 01/2013 : Adapted for large output files with at least 2 variables > 2 GiB
	17
	18	implicit none
	19
	20	include "netcdf.inc" ! NetCDF definitions
	21
	22	character(len=128) :: infile ! input file name (diagfi.nc or stats.nc format)
	23	character(len=128) :: outfile ! output file name
	24	character(len=64) :: text ! to store some text
	25	character(len=64) :: tmpvarname ! temporarily store a variable name
	26	integer :: infid ! NetCDF input file ID
	27	integer :: outfid ! NetCDF output file ID
	28	integer :: tmpvarid ! temporarily store a variable ID
	29	integer :: tmpdimid ! temporarily store a dimension ID
	30	integer :: tmpndims ! temporarily store # of dimensions of a variable
	31	integer :: nbvarinfile ! # of variables in input file
	32	integer :: nbattr ! # of attributes of a given variable in input file
	33	integer :: nbvar3dinfile ! # of 3D (lon,lat,time) variables in input file
	34	integer :: nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file
	35
	36	integer :: i,j
	37	integer :: lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
	38	integer :: lon_varid,lat_varid,alt_varid,time_varid ! NetCDF variable IDs
	39	!integer gcm_layers_dimid ! NetCDF dimension ID for # of layers in GCM
	40	integer :: sigma_varid,aps_varid,bps_varid
	41	integer :: phisinit_varid
	42	integer,dimension(4) :: datashape ! shape of 4D datasets
	43	integer :: ierr ! NetCDF routines return code
	44	character (len=64), dimension(:), allocatable :: var
	45	! var(): names of variables that will be processed
	46	integer nbvar ! # of variables to process
	47	integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile)
	48	real :: miss_val=-9.99e+33 ! special "missing value" to specify missing data
	49	real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
	50	real,dimension(:),allocatable :: inlat ! input latitude
	51	integer :: inlatlength ! # of elements in input latitude
	52	real,dimension(:),allocatable :: inlon ! input longitude
	53	integer :: inlonlength ! # of elements in input longitude
	54	real,dimension(:),allocatable :: alt ! altitude
	55	integer :: altlength ! # of elements in altitude
	56	real,dimension(:),allocatable :: time ! input time
	57	integer :: timelength ! # of elements in time(:)
	58	real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
	59	real,dimension(:),allocatable :: sigma ! sigma levels
	60	real,dimension(:,:),allocatable :: inphisinit ! input ground geopotential
	61	real,dimension(:),allocatable :: lon ! output longitude
	62	integer :: lonlength ! # of elements in lon()
	63	real,dimension(:),allocatable :: wklon ! work longitude (with modulo element)
	64	integer :: wklonlength ! # of elements in wklon()
	65	real,dimension(:),allocatable :: lat ! output latitude
	66	integer :: latlength ! # of elements in lat()
	67	real,dimension(:),allocatable :: wklat ! work latitude (includes poles)
	68	integer :: wklatlength ! # of elements in wklat()
	69	!real,dimension(:),allocatable :: lon_bound ! output longitude boundaries
	70	!real,dimension(:),allocatable :: lat_bound ! output latitude boundaries
	71	real,dimension(:),allocatable :: in_lon_bound ! input longitude boundaries
	72	real,dimension(:),allocatable :: in_lat_bound ! input latitude boundaries
	73	real,dimension(:),allocatable :: wk_lon_bound ! work longitude boundaries
	74	real,dimension(:),allocatable :: wk_lat_bound ! work latitude boundaries
	75	real,dimension(:,:),allocatable :: in_2d_data ! input 2D (lon-lat) dataset
	76	real,dimension(:,:),allocatable :: wk_2d_data ! work 2D dataset
	77	real,dimension(:,:),allocatable :: out_2d_data ! output 2D dataset
	78	real,dimension(:,:,:),allocatable :: in_3d_data ! intput 3D dataset
	79	real,dimension(:,:,:),allocatable :: wk_3d_data ! work 3D dataset
	80	real,dimension(:,:,:),allocatable :: out_3d_data ! output 3D dataset
	81	real,dimension(:,:,:,:),allocatable :: in_4d_data ! intput 4D dataset
	82	real,dimension(:,:,:,:),allocatable :: wk_4D_data ! work 4D dataset
	83	real,dimension(:,:,:,:),allocatable :: out_4d_data ! output 4D dataset
	84
	85	real :: pi ! =3.14...
	86	logical :: have_sigma ! Flag: true if sigma levels are known (false if hybrid
	87	! coordinates are used)
	88	logical :: have_geopot ! Flag: true if input file contains ground geopotential
	89	! phisinit()
	90	logical :: out_mod_lon ! Flag: true if output grid has modulo longitude (ie:
	91	! first and last point are in fact at same longitude)
	92	logical :: out_has_poles ! Flag: true if output grid includes North and South
	93	! poles
	94
	95	integer, dimension(4) :: edges,corner ! needed to write variables for big files
	96
	97	!===============================================================================
	98	! 1. Input parameters
	99	!===============================================================================
	100	pi=2.*asin(1.)
	101
	102	!===============================================================================
	103	! 1.1 Input file
	104	!===============================================================================
	105
	106	write(,) ""
	107	write(,) " Program valid for diagfi.nc, concatnc.nc and stats.nc files"
	108	write(,) "Enter input file name:"
	109
	110	read(*,'(a128)') infile
	111	write(,) ""
	112
	113	! open input file
	114
	115	ierr = NF_OPEN(infile,NF_NOWRITE,infid)
	116	if (ierr.ne.NF_NOERR) then
	117	write(,) 'ERROR: Pb opening file ',trim(infile)
	118	stop ""
	119	endif
	120
	121	!===============================================================================
	122	! 1.2 Get # and names of variables in input file
	123	!===============================================================================
	124
	125	ierr=NF_INQ_NVARS(infid,nbvarinfile)
	126	if (ierr.ne.NF_NOERR) then
	127	write(,) 'ERROR: Failed geting number of variables from file'
	128	stop
	129	endif
	130
	131	write(,)" The following variables have been found:"
	132	nbvar3dinfile=0
	133	nbvar4dinfile=0
	134	do i=1,nbvarinfile
	135	! get name of variable # i
	136	ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
	137	! check if it is a 3D variable
	138	ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
	139	if (tmpndims.eq.3) then
	140	nbvar3dinfile=nbvar3dinfile+1
	141	write(,) trim(tmpvarname)
	142	endif
	143	! check if it is a 4D variable
	144	ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
	145	if (tmpndims.eq.4) then
	146	nbvar4dinfile=nbvar4dinfile+1
	147	write(,) trim(tmpvarname)
	148	endif
	149	enddo
	150
	151	allocate(var(nbvar3dinfile+nbvar4dinfile))
	152
	153	write(,) ""
	154	write(,) "Which variable do you want to keep?"
	155	write(,) "all or list of <variables> (separated by <Return>s)"
	156	write(,) "(an empty line , i.e: just <Return>, implies end of list)"
	157	nbvar=0
	158	read(*,'(a64)') tmpvarname
	159	do while ((tmpvarname.ne.' ').and.(trim(tmpvarname).ne.'all'))
	160	! check if tmpvarname is valid
	161	ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
	162	if (ierr.eq.NF_NOERR) then ! valid name
	163	nbvar=nbvar+1
	164	var(nbvar)=tmpvarname
	165	else ! invalid name
	166	write(,) 'Error: ',trim(tmpvarname),' is not a valid name'
	167	write(,) ' (we''ll skip that one)'
	168	endif
	169	read(*,'(a64)') tmpvarname
	170	enddo
	171
	172	! handle "all" case
	173	if (tmpvarname.eq.'all') then
	174	nbvar=0
	175	do i=1,nbvarinfile
	176	! look for 4D variables
	177	ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
	178	if (tmpndims.eq.4) then
	179	nbvar=nbvar+1
	180	! get the corresponding name
	181	ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
	182	var(nbvar)=tmpvarname
	183	endif
	184	enddo
	185	endif
	186
	187	! Check that there is at least 1 variable to process
	188	if (nbvar.eq.0) then
	189	write(,) 'No variables to process !?'
	190	write(,) 'Might as well stop here'
	191	stop ""
	192	else
	193	write(,) ""
	194	write(,) 'OK, the following variables will be processed:'
	195	do i=1,nbvar
	196	write(,) var(i)
	197	enddo
	198	endif
	199
	200	!===============================================================================
	201	! 1.3 Get input grids in lon,lat,alt,time,
	202	! as well as hybrid coordinates aps() and bps() (or sigma levels sigma())
	203	! and eventually phisinit() from input file
	204	!===============================================================================
	205
	206	! 1.3.1 input longitude, latitude, altitude and time
	207
	208	! latitude
	209	ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid)
	210	if (ierr.ne.NF_NOERR) then
	211	stop "Error: Failed to get latitude dimension ID"
	212	else
	213	ierr=NF_INQ_VARID(infid,"latitude",tmpvarid)
	214	if (ierr.ne.NF_NOERR) then
	215	stop "Error: Failed to get latitude ID"
	216	else
	217	ierr=NF_INQ_DIMLEN(infid,tmpdimid,inlatlength)
	218	if (ierr.ne.NF_NOERR) then
	219	stop "Error: Failed to get latitude length"
	220	else
	221	allocate(inlat(inlatlength))
	222	ierr=NF_GET_VAR_REAL(infid,tmpvarid,inlat)
	223	if (ierr.ne.NF_NOERR) then
	224	stop "Error: Failed reading latitude"
	225	endif
	226	endif
	227	endif
	228	endif
	229
	230	! check that these latitudes are 'GCM' latitudes (i.e. poles are included)
	231	if ((abs(inlat(1)-90.0)).gt.0.001) then
	232	write(,) "Error: Input latitudes should include north pole, but"
	233	write(,) " lat(1)=",inlat(1)
	234	stop
	235	endif
	236	if ((abs(inlat(inlatlength)+90.0)).gt.0.001) then
	237	write(,) "Error: Input latitudes should include south pole, but"
	238	write(,) " lat(inlatlength)=",inlat(inlatlength)
	239	stop
	240	endif
	241
	242
	243	! longitude
	244	ierr=NF_INQ_DIMID(infid,"longitude",tmpdimid)
	245	if (ierr.ne.NF_NOERR) then
	246	stop "Error: Failed to get longitude dimension ID"
	247	else
	248	ierr=NF_INQ_VARID(infid,"longitude",tmpvarid)
	249	if (ierr.ne.NF_NOERR) then
	250	stop "Error: Failed to get longitude ID"
	251	else
	252	ierr=NF_INQ_DIMLEN(infid,tmpdimid,inlonlength)
	253	if (ierr.ne.NF_NOERR) then
	254	stop "Error: Failed to get longitude length"
	255	else
	256	allocate(inlon(inlonlength))
	257	ierr=NF_GET_VAR_REAL(infid,tmpvarid,inlon)
	258	if (ierr.ne.NF_NOERR) then
	259	stop "Error: Failed reading longitude"
	260	endif
	261	endif
	262	endif
	263	endif
	264
	265	! check that these longitudes are 'GCM' longitudes (i.e range from -180 to 180)
	266	if (abs(180.+inlon(1)).gt.0.001) then
	267	write(,) "Error: Input latitudes should start at -180, but"
	268	write(,) " lon(1)=",inlon(1)
	269	stop
	270	endif
	271	if (abs(inlon(inlonlength)-180).gt.0.001) then
	272	write(,) "Error: Input latitudes should end at 180, but"
	273	write(,) " lon(inlonlength)=",inlon(inlonlength)
	274	stop
	275	endif
	276
	277
	278	! altitude
	279	ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid)
	280	if (ierr.ne.NF_NOERR) then
	281	stop "Error: Failed to get altitude dimension ID"
	282	else
	283	ierr=NF_INQ_VARID(infid,"altitude",tmpvarid)
	284	if (ierr.ne.NF_NOERR) then
	285	stop "Error: Failed to get altitude ID"
	286	else
	287	ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
	288	if (ierr.ne.NF_NOERR) then
	289	stop "Error: Failed to get altitude length"
	290	else
	291	allocate(alt(altlength))
	292	ierr=NF_GET_VAR_REAL(infid,tmpvarid,alt)
	293	if (ierr.ne.NF_NOERR) then
	294	stop "Error: Failed reading altitude"
	295	endif
	296	endif
	297	endif
	298	endif
	299
	300	! time
	301	ierr=NF_INQ_DIMID(infid,"Time",tmpdimid)
	302	if (ierr.ne.NF_NOERR) then
	303	stop "Error: Failed to get Time dimension ID"
	304	else
	305	ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
	306	if (ierr.ne.NF_NOERR) then
	307	stop "Error: Failed to get Time ID"
	308	else
	309	ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
	310	if (ierr.ne.NF_NOERR) then
	311	stop "Error: Failed to get Time length"
	312	else
	313	allocate(time(timelength))
	314	ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
	315	if (ierr.ne.NF_NOERR) then
	316	stop "Error: Failed reading Time"
	317	endif
	318	endif
	319	endif
	320	endif
	321
	322	! 1.3.2 Get hybrid coordinates (or sigma levels)
	323
	324	! start by looking for sigma levels
	325	ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
	326	if (ierr.ne.NF_NOERR) then
	327	have_sigma=.false.
	328	write(,) "Could not find sigma levels... will look for hybrid coordinates"
	329	else
	330	have_sigma=.true.
	331	allocate(sigma(altlength))
	332	ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
	333	if (ierr.ne.NF_NOERR) then
	334	stop "Error: Failed reading sigma"
	335	endif
	336	endif
	337
	338	! if no sigma levels, look for hybrid coordinates
	339	if (.not.have_sigma) then
	340	! hybrid coordinate aps
	341	ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
	342	if (ierr.ne.NF_NOERR) then
	343	stop "Error: Failed to get aps ID"
	344	else
	345	allocate(aps(altlength))
	346	ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
	347	if (ierr.ne.NF_NOERR) then
	348	stop "Error: Failed reading aps"
	349	endif
	350	endif
	351
	352	! hybrid coordinate bps
	353	ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
	354	if (ierr.ne.NF_NOERR) then
	355	stop "Error: Failed to get bps ID"
	356	else
	357	allocate(bps(altlength))
	358	ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
	359	if (ierr.ne.NF_NOERR) then
	360	stop "Error: Failed reading bps"
	361	endif
	362	endif
	363	endif !of if (.not.have_sigma)
	364
	365	! 1.3.3 Get ground geopotential phisinit, if available
	366
	367	! look for 'phisinit' in current file
	368	ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
	369	if (ierr.ne.NF_NOERR) then
	370	write(,) "Warning: Failed to get phisinit ID from file ",trim(infile)
	371	write(,) " ...will not store geopotential in output... "
	372	have_geopot=.false.
	373	else
	374	have_geopot=.true.
	375	! Get input physinit
	376	allocate(inphisinit(inlonlength,inlatlength))
	377	ierr=NF_GET_VAR_REAL(infid,tmpvarid,inphisinit)
	378	if (ierr.ne.NF_NOERR) then
	379	stop "Error: Failed reading phisinit"
	380	endif
	381	endif
	382
	383	! 1.3.4 Create input longitude and latitude boundaries
	384
	385	! build input longitude boundaries (in radians)
	386	allocate(in_lon_bound(inlonlength))
	387	do i=1,inlonlength-1
	388	in_lon_bound(i)=0.5(inlon(i+1)+inlon(i))pi/180.0
	389	enddo
	390	! we have inlon(1)=inlon(inlonlength) modulo 360
	391	! and thus in_lon_bound(inlonlength)=in_lon_bound(1)+360
	392	in_lon_bound(inlonlength)=2.*pi+in_lon_bound(1)
	393	!do i=1,inlonlength
	394	! write(,) "i=",i,"180/piin_lon_bound(i)=",(180./pi)in_lon_bound(i)
	395	!enddo
	396
	397	! build input latitude boundaries (in radians)
	398	allocate(in_lat_bound(inlatlength-1))
	399	do i=1,inlatlength-1
	400	in_lat_bound(i)=0.5(inlat(i)+inlat(i+1))pi/180.0
	401	enddo
	402	!do i=1,inlatlength-1
	403	! write(,) "i=",i,"180/piin_lat_bound(i)",(180./pi)in_lat_bound(i)
	404	!enddo
	405
	406	!===============================================================================
	407	! 1.4 Get output longitude and latitude coordinates
	408	!===============================================================================
	409
	410	write(,) ""
	411	write(,) "Output horizontal grid:"
	412	write(,) "Number of grid points in longitude?"
	413	read(,) lonlength
	414	write(,) "Enter longitudes (degrees, in [-180:180]), "
	415	write(,) " in increasing order (one per line):"
	416	allocate(lon(lonlength))
	417	do i=1,lonlength
	418	read(,) lon(i)
	419	enddo
	420
	421	!! build 'work' longitude (which must be a modulo axis; i.e. first and
	422	!! last points e.g. -180 and 180 are included)
	423	if (abs((lon(1)+360.-lon(lonlength))).le.0.01) then
	424	! the axis already has modulo endpoints
	425	out_mod_lon=.true.
	426	wklonlength=lonlength
	427	allocate(wklon(wklonlength))
	428	wklon(1:lonlength)=lon(1:lonlength)
	429	else
	430	! add an extra point
	431	out_mod_lon=.false.
	432	wklonlength=lonlength+1
	433	allocate(wklon(wklonlength))
	434	wklon(1:lonlength)=lon(1:lonlength)
	435	wklon(wklonlength)=wklon(1)+360.0
	436	endif
	437
	438	! build work longitude boundaries (in radians)
	439	allocate(wk_lon_bound(wklonlength))
	440	do i=1,lonlength-1
	441	wk_lon_bound(i)=0.5(wklon(i+1)+wklon(i))pi/180.0
	442	enddo
	443	if (out_mod_lon) then
	444	! we have lon(1)=lon(lonlength) modulo 360
	445	! and thus lon_bound(lonlength)=lon_bound(1)+360
	446	wk_lon_bound(wklonlength)=2.*pi+wk_lon_bound(1)
	447	else
	448	wk_lon_bound(wklonlength-1)=0.5(wklon(wklonlength)+wklon(wklonlength-1))pi/180.0
	449	wk_lon_bound(wklonlength)=2.*pi+wk_lon_bound(1)
	450	endif
	451
	452
	453	write(,) "Number of grid points in latitude?"
	454	read(,) latlength
	455	write(,) "Enter latitudes (degrees), in decreasing order (from northernmost"
	456	write(,) " to southernmost), one per line:"
	457	allocate(lat(latlength))
	458	do i=1,latlength
	459	read(,) lat(i)
	460	enddo
	461
	462	! build 'work' latitude (which must include poles, ie lat=90 and -90)
	463	if (abs(lat(1)-90.0).le.0.001) then
	464	out_has_poles=.true.
	465	wklatlength=latlength
	466	allocate(wklat(wklatlength))
	467	wklat(1:latlength)=lat(1:latlength)
	468	else
	469	out_has_poles=.false.
	470	! add poles
	471	wklatlength=latlength+2
	472	allocate(wklat(wklatlength))
	473	wklat(1)=90
	474	wklat(2:latlength+1)=lat(1:latlength)
	475	wklat(wklatlength)=-90
	476	endif
	477
	478	! build work latitude boundaries (in radians)
	479	allocate(wk_lat_bound(wklatlength-1))
	480	if (out_has_poles) then
	481	do i=1,wklatlength-1
	482	wk_lat_bound(i)=0.5(wklat(i)+wklat(i+1))pi/180.0
	483	enddo
	484	else
	485	! put northermost boundary near pole
	486	wk_lat_bound(1)=(90-0.01(90.-lat(1)))pi/180.0
	487	do i=2,wklatlength-2
	488	wk_lat_bound(i)=0.5(wklat(i)+wklat(i+1))pi/180.0
	489	enddo
	490	! put southernmost boundary near pole
	491	wk_lat_bound(wklatlength-1)=(-90.0-0.01(-90.-lat(latlength)))pi/180.0
	492	endif
	493
	494	!do i=1,wklatlength-1
	495	! write(,) "i=",i,"180/piwk_lat_bound(i)",(180./pi)wk_lat_bound(i)
	496	!enddo
	497
	498	!===============================================================================
	499	! 1.5 Output file
	500	!===============================================================================
	501	write(,) ""
	502	outfile=infile(1:len_trim(infile)-3)//"_h.nc"
	503	write(,) "Output file name is: ",trim(outfile)
	504
	505
	506	!===============================================================================
	507	! 2. Create output file and initialize definitions of variables and dimensions
	508	!===============================================================================
	509
	510	!===============================================================================
	511	! 2.1. Output file
	512	!===============================================================================
	513
	514	! Create output file
	515	ierr=NF_CREATE(outfile,IOR(NF_CLOBBER,NF_64BIT_OFFSET),outfid)
	516	if (ierr.ne.NF_NOERR) then
	517	write(,)"Error: could not create file ",outfile
	518	stop
	519	endif
	520
	521	!===============================================================================
	522	! 2.2. Define dimensions
	523	!===============================================================================
	524	! longitude
	525	ierr=NF_DEF_DIM(outfid,"longitude",lonlength,lon_dimid)
	526	if (ierr.ne.NF_NOERR) then
	527	stop "Error: Could not define longitude dimension"
	528	endif
	529
	530	! latitude
	531	ierr=NF_DEF_DIM(outfid,"latitude",latlength,lat_dimid)
	532	if (ierr.ne.NF_NOERR) then
	533	stop "Error: Could not define latitude dimension"
	534	endif
	535
	536	! altitude
	537	ierr=NF_DEF_DIM(outfid,"altitude",altlength,alt_dimid)
	538	if (ierr.ne.NF_NOERR) then
	539	stop "Error: Could not define altitude dimension"
	540	endif
	541
	542	! time
	543	ierr=NF_DEF_DIM(outfid,"Time",NF_UNLIMITED,time_dimid)
	544	if (ierr.ne.NF_NOERR) then
	545	stop "Error: Could not define latitude dimension"
	546	endif
	547
	548	!===============================================================================
	549	! 2.3. Define variables and their attributes
	550	!===============================================================================
	551
	552	! 2.3.1 Define 1D variables
	553
	554	! longitude
	555	datashape(1)=lon_dimid
	556	ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid)
	557	if (ierr.ne.NF_NOERR) then
	558	stop "Error: Could not define longitude variable"
	559	endif
	560
	561	! longitude attributes
	562	text='east longitude'
	563	ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'long_name',len_trim(text),text)
	564	if (ierr.ne.NF_NOERR) then
	565	stop "Error: Problem writing long_name for longitude"
	566	endif
	567	text='degrees_east'
	568	ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'units',len_trim(text),text)
	569	if (ierr.ne.NF_NOERR) then
	570	stop "Error: Problem writing units for longitude"
	571	endif
	572
	573	! latitude
	574	datashape(2)=lat_dimid
	575	ierr=NF_DEF_VAR(outfid,"latitude",NF_REAL,1,datashape(2),lat_varid)
	576	if (ierr.ne.NF_NOERR) then
	577	stop "Error: Could not define latitude variable"
	578	endif
	579
	580	! latitude attributes
	581	text='north latitude'
	582	ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'long_name',len_trim(text),text)
	583	if (ierr.ne.NF_NOERR) then
	584	stop "Error: Problem writing long_name for latitude"
	585	endif
	586	text='degrees_north'
	587	ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'units',len_trim(text),text)
	588	if (ierr.ne.NF_NOERR) then
	589	stop "Error: Problem writing units for latitude"
	590	endif
	591
	592	! altitude
	593	datashape(3)=alt_dimid
	594	ierr=NF_DEF_VAR(outfid,"altitude",NF_REAL,1,datashape(3),alt_varid)
	595	if (ierr.ne.NF_NOERR) then
	596	stop "Error: Could not define altitude variable"
	597	endif
	598
	599	! altitude attributes
	600	! preliminary stuff, get input file altitude variable ID
	601	ierr=NF_INQ_VARID(infid,"altitude",tmpvarid)
	602
	603	! look for a 'long_name' attribute
	604	text=" "
	605	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,"long_name",text)
	606	if (ierr.eq.NF_NOERR) then
	607	! found the attribute; write it to output file
	608	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
	609	endif
	610
	611	! look for a 'unit' attribute
	612	text=" "
	613	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,"units",text)
	614	if (ierr.eq.NF_NOERR) then
	615	! found the attribute; write it to output file
	616	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
	617	endif
	618
	619	! look for a 'positive' attribute
	620	text=" "
	621	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,"positive",text)
	622	if (ierr.eq.NF_NOERR) then
	623	! found the attribute; write it to output file
	624	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
	625	endif
	626
	627	! sigma levels or hybrid coordinates
	628	if (have_sigma) then
	629	ierr=NF_DEF_VAR(outfid,"sigma",NF_REAL,1,alt_dimid,sigma_varid)
	630	if (ierr.ne.NF_NOERR) then
	631	stop "Error: Could not define sigma variable"
	632	endif
	633	else ! hybrid coordinates
	634	ierr=NF_DEF_VAR(outfid,"aps",NF_REAL,1,alt_dimid,aps_varid)
	635	if (ierr.ne.NF_NOERR) then
	636	stop "Error: Could not define aps variable"
	637	endif
	638	ierr=NF_DEF_VAR(outfid,"bps",NF_REAL,1,alt_dimid,bps_varid)
	639	if (ierr.ne.NF_NOERR) then
	640	stop "Error: Could not define bps variable"
	641	endif
	642	endif
	643
	644	! sigma levels (or hybrid coordinates) attributes
	645	if (have_sigma) then
	646	text="sigma levels"
	647	ierr=NF_PUT_ATT_TEXT(outfid,sigma_varid,'long_name',len_trim(text),text)
	648	if (ierr.ne.NF_NOERR) then
	649	stop "Error: Problem writing long_name for sigma"
	650	endif
	651	else ! hybrid coordinates
	652	text="hybrid pressure at midlayers"
	653	ierr=NF_PUT_ATT_TEXT(outfid,aps_varid,'long_name',len_trim(text),text)
	654	if (ierr.ne.NF_NOERR) then
	655	stop "Error: Problem writing long_name for aps"
	656	endif
	657	text="hybrid sigma at midlayers"
	658	ierr=NF_PUT_ATT_TEXT(outfid,bps_varid,'long_name',len_trim(text),text)
	659	if (ierr.ne.NF_NOERR) then
	660	stop "Error: Problem writing long_name for bps"
	661	endif
	662	endif ! of if (have_sigma)
	663
	664	! time
	665	datashape(4)=time_dimid
	666	ierr=NF_DEF_VAR(outfid,"Time",NF_REAL,1,datashape(4),time_varid)
	667	if (ierr.ne.NF_NOERR) then
	668	stop "Error: Could not define Time variable"
	669	endif
	670
	671	! time attributes
	672	! preliminary stuff, get input file time variable ID
	673	ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
	674	! look for a 'unit' attribute
	675	text=" "
	676	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,"units",text)
	677	if (ierr.eq.NF_NOERR) then
	678	! found the attribute; write it to output file
	679	ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
	680	else
	681	! write something not too stupid
	682	text='days since 0000-01-1 00:00:00'
	683	ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
	684	if (ierr.ne.NF_NOERR) then
	685	stop "Error: Problem writing units for Time"
	686	endif
	687	endif
	688
	689	! 2.3.2 Define 2D (lon-lat, time independent) variables
	690
	691	! ground geopotential
	692	if (have_geopot) then
	693	if (have_geopot) then
	694	ierr=NF_DEF_VAR(outfid,"phisinit",NF_REAL,2,datashape,phisinit_varid)
	695	if (ierr.ne.NF_NOERR) then
	696	stop "Error: Could not define phisinit variable"
	697	endif
	698	endif
	699	! ground geopotential attributes
	700	text='Geopotential at the surface'
	701	ierr=NF_PUT_ATT_TEXT(outfid,phisinit_varid,'long_name',len_trim(text),text)
	702	if (ierr.ne.NF_NOERR) then
	703	stop "Error: Problem writing long_name for phisinit"
	704	endif
	705	endif ! of if (have_geopot)
	706
	707	! 2.3.3 Define 3D&4D variables (ie: surface/volume+time variables)
	708
	709	! variables requested by user
	710	allocate(var_id(nbvar))
	711	do i=1,nbvar
	712	! Get the (input file) ID of the variable (stored in tmpvarid)
	713	ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
	714	if (ierr.ne.NF_NOERR) then
	715	write(,) 'Error, failed to get ID for input variable ',trim(var(i))
	716	stop ""
	717	endif
	718
	719	! Get # of dimensions of the variable
	720	! and define the variable in output file
	721	write(,) ""
	722	write(,) "Creating ",trim(var(i))
	723	ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
	724	if (tmpndims.eq.3) then
	725	datashape(3)=time_dimid
	726	ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,3,datashape,var_id(i))
	727	else
	728	if (tmpndims.eq.4) then
	729	datashape(3)=alt_dimid
	730	datashape(4)=time_dimid
	731	ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
	732	else
	733	write(,) "Error: number of dimensions of input variable ",trim(var(i))
	734	write(,) " is ",tmpndims
	735	stop
	736	endif
	737	endif
	738	if (ierr.ne.NF_NOERR) then
	739	write(,) 'Error, could not define variable ',trim(var(i))
	740	stop ""
	741	endif
	742
	743	! Get the (input file) ID for the variable and
	744	! the # of attributes associated to that variable
	745	ierr=NF_INQ_VARNATTS(infid,tmpvarid,nbattr)
	746	if (ierr.ne.NF_NOERR) then
	747	write(,) 'Error, could not get number of attributes for variable ',&
	748	trim(var(i))
	749	stop ""
	750	endif
	751	! inititialize j == number of attributes written to output
	752	j=0
	753
	754	! look for a "long_name" attribute (or eventually a 'title' attribute)
	755	text=' '
	756	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'long_name',text)
	757	if (ierr.ne.NF_NOERR) then ! no long_name attribute
	758	! try to find an equivalent 'title' attribute
	759	text=' '
	760	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'title',text)
	761	if (ierr.eq.NF_NOERR) then ! found 'title' attribute
	762	write(,) "Found title ",trim(text)
	763	j=j+1
	764	! write it as a 'long_name' attribute
	765	ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
	766	if (ierr.ne.NF_NOERR) then
	767	write(,) "Error failed to copy title attribute:",trim(text)
	768	stop ""
	769	endif
	770	endif
	771	else ! found long_name; write it to outfile
	772	write(,) "Found long_name ",trim(text)
	773	j=j+1
	774	ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
	775	if (ierr.ne.NF_NOERR) then
	776	write(,) "Error failed to copy long_name attribute:",trim(text)
	777	stop""
	778	endif
	779	endif
	780
	781	! look for a "units" attribute
	782	text=' '
	783	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',text)
	784	if (ierr.eq.NF_NOERR) then ! found 'units' attribute
	785	write(,) "Found units ",trim(text)
	786	j=j+1
	787	! write it to output
	788	ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'units',len_trim(text),text)
	789	if (ierr.ne.NF_NOERR) then
	790	write(,) "Error failed to copy units attribute:",trim(text)
	791	stop""
	792	endif
	793	endif
	794
	795	! look for a "missing_value" attribute
	796	ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val)
	797	if (ierr.eq.NF_NOERR) then ! found 'missing_value' attribute
	798	write(,) "Found missing_value ",miss_val
	799	j=j+1
	800	else ! no 'missing_value' attribute, set miss_val to default
	801	miss_val=miss_val_def
	802	endif
	803
	804	! write the missing_value attribute to output
	805	ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val)
	806	if (ierr.ne.NF_NOERR) then
	807	stop "Error, failed to write missing_value attribute"
	808	endif
	809
	810	! warn if some attributes were missed
	811	if (j.ne.nbattr) then
	812	write(,)'Warning, it seems some attributes of variable ',trim(var(i))
	813	write(,)"were not transfered to the new file"
	814	write(,)'nbattr:',nbattr,' j:',j
	815	endif
	816
	817	enddo ! of do i=1,nbvar
	818
	819	!===============================================================================
	820	! 2.4. Write dimensions (and time-independent variables)
	821	!===============================================================================
	822	! Switch out of NetCDF define mode
	823	ierr=NF_ENDDEF(outfid)
	824	if (ierr.ne.NF_NOERR) then
	825	stop "Error: Could not switch out of define mode"
	826	endif
	827
	828	! Write longitude
	829	ierr=NF_PUT_VAR_REAL(outfid,lon_varid,lon)
	830	if (ierr.ne.NF_NOERR) then
	831	stop "Error: Could not write longitude data to output file"
	832	endif
	833
	834	! Write latitude
	835	ierr=NF_PUT_VAR_REAL(outfid,lat_varid,lat)
	836	if (ierr.ne.NF_NOERR) then
	837	stop "Error: Could not write latitude data to output file"
	838	endif
	839
	840	! write altitude
	841	ierr=NF_PUT_VAR_REAL(outfid,alt_varid,alt)
	842	if (ierr.ne.NF_NOERR) then
	843	stop "Error: Could not write altitude data to output file"
	844	endif
	845
	846	! Write sigma levels (or hybrid coordinates)
	847	if (have_sigma) then
	848	ierr=NF_PUT_VAR_REAL(outfid,sigma_varid,sigma)
	849	if (ierr.ne.NF_NOERR) then
	850	stop "Error: Could not write sigma data to output file"
	851	endif
	852	else ! hybrid coordinates
	853	ierr=NF_PUT_VAR_REAL(outfid,aps_varid,aps)
	854	if (ierr.ne.NF_NOERR) then
	855	stop "Error: Could not write aps data to output file"
	856	endif
	857	ierr=NF_PUT_VAR_REAL(outfid,bps_varid,bps)
	858	if (ierr.ne.NF_NOERR) then
	859	stop "Error: Could not write bps data to output file"
	860	endif
	861	endif
	862
	863	! write time
	864	ierr=NF_PUT_VARA_REAL(outfid,time_varid,1,timelength,time)
	865	if (ierr.ne.NF_NOERR) then
	866	stop "Error: Could not write Time data to output file"
	867	endif
	868
	869	!===============================================================================
	870	! 3. Interpolate and write 2D variables
	871	!===============================================================================
	872	write(,) "interpolate 2D variables"
	873	allocate(in_2d_data(inlonlength,inlatlength))
	874	allocate(wk_2d_data(wklonlength,wklatlength))
	875	allocate(out_2d_data(lonlength,latlength))
	876
	877	! ground geopotential
	878	if (have_geopot) then
	879	! load input geopotential: get ID for input phisinit
	880	ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
	881	if (ierr.ne.NF_NOERR) then
	882	stop "Error: Failed to get phisinit ID"
	883	endif
	884	! Get physinit
	885	ierr=NF_GET_VAR_REAL(infid,tmpvarid,in_2d_data)
	886	if (ierr.ne.NF_NOERR) then
	887	stop "Error: Failed reading input phisinit"
	888	endif
	889
	890	! interpolate onto new grid
	891	call interp_horiz(in_2d_data,wk_2d_data,inlonlength-1,inlatlength-1,&
	892	wklonlength-1,wklatlength-1,1,&
	893	in_lon_bound,in_lat_bound,wk_lon_bound,wk_lat_bound)
	894
	895	! copy (and possibly reshape) data to output grid
	896	if (out_has_poles) then
	897	out_2d_data(1:lonlength,1:latlength)=wk_2d_data(1:lonlength,1:latlength)
	898	else
	899	out_2d_data(1:lonlength,1:latlength)=wk_2d_data(1:lonlength,2:latlength+1)
	900	endif
	901
	902	! write interpolated phisinit to output
	903	ierr=NF_PUT_VAR_REAL(outfid,phisinit_varid,out_2d_data)
	904	if (ierr.ne.NF_NOERR) then
	905	stop "Error: Could not write phisinit data to output file"
	906	endif
	907	endif ! of if (have_geopot)
	908
	909	!===============================================================================
	910	! 4. Interpolate and write 3D (surface+time) variables
	911	!===============================================================================
	912	write(,) "interpolate 3D variables"
	913	allocate(in_3d_data(inlonlength,inlatlength,timelength))
	914	allocate(wk_3d_data(wklonlength,wklatlength,timelength))
	915	allocate(out_3d_data(lonlength,latlength,timelength))
	916
	917	do i=1,nbvar ! loop on all selected 3D&4D variables
	918	! get input file ID for this variable
	919	ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
	920	if (ierr.ne.NF_NOERR) then
	921	write(,) "Error: Failed to get input ID for ",trim(var(i))
	922	stop
	923	endif
	924	! get # of dimensions for this variable
	925	ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
	926
	927	if (tmpndims.eq.3) then ! if it is indeed at 3D variable
	928	! get the data
	929	ierr=NF_GET_VAR_REAL(infid,tmpvarid,in_3d_data)
	930	if (ierr.ne.NF_NOERR) then
	931	write(,) "Error: Failed reading input ",trim(var(i))
	932	stop
	933	endif
	934	! interpolate data
	935	do j=1,timelength
	936	call interp_horiz(in_3d_data(1,1,j),wk_3d_data(1,1,j),&
	937	inlonlength-1,inlatlength-1,&
	938	wklonlength-1,wklatlength-1,1,&
	939	in_lon_bound,in_lat_bound,wk_lon_bound,wk_lat_bound)
	940	! copy (and possibly reshape) data to output grid
	941	if (out_has_poles) then
	942	out_3d_data(1:lonlength,1:latlength,j)=&
	943	wk_3d_data(1:lonlength,1:latlength,j)
	944	else
	945	out_3d_data(1:lonlength,1:latlength,j)=&
	946	wk_3d_data(1:lonlength,2:latlength+1,j)
	947	endif
	948	enddo
	949	! write interpolated data to output
	950	corner(:)=1
	951	edges(1)=lonlength
	952	edges(2)=latlength
	953	edges(3)=timelength
	954	ierr=NF_PUT_VARA_REAL(outfid,var_id(i),corner(1:3),edges(1:3),out_3d_data)
	955	if (ierr.ne.NF_NOERR) then
	956	write(,) "Error: Could not write ",trim(var(i))," data to output file"
	957	stop
	958	else
	959	write(,) "wrote ",trim(var(i))
	960	endif
	961	endif ! of if (tmpndims.eq.3)
	962
	963	enddo ! of do i=1,nbvar
	964
	965	!===============================================================================
	966	! 5. Interpolate and write 4D variables
	967	!===============================================================================
	968	allocate(in_4d_data(inlonlength,inlatlength,altlength,timelength))
	969	allocate(wk_4d_data(wklonlength,wklatlength,altlength,timelength))
	970	allocate(out_4d_data(lonlength,latlength,altlength,timelength))
	971
	972	do i=1,nbvar ! loop on all selected 3D&4D variables
	973	! get input file ID for this variable
	974	ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
	975	if (ierr.ne.NF_NOERR) then
	976	write(,) "Error: Failed to get input ID for ",trim(var(i))
	977	stop
	978	endif
	979	! get # of dimensions for this variable
	980	ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
	981
	982	if (tmpndims.eq.4) then ! if it is indeed at 4D variable
	983	! get the data
	984	ierr=NF_GET_VAR_REAL(infid,tmpvarid,in_4d_data)
	985	if (ierr.ne.NF_NOERR) then
	986	write(,) "Error: Failed reading input ",trim(var(i))
	987	stop
	988	endif
	989	! interpolate data
	990	do j=1,timelength
	991	call interp_horiz(in_4d_data(1,1,1,j),wk_4d_data(1,1,1,j),&
	992	inlonlength-1,inlatlength-1,&
	993	wklonlength-1,wklatlength-1,altlength,&
	994	in_lon_bound,in_lat_bound,wk_lon_bound,wk_lat_bound)
	995	! copy (and possibly reshape) data to output grid
	996	if (out_has_poles) then
	997	out_4d_data(1:lonlength,1:latlength,1:altlength,j)=&
	998	wk_4d_data(1:lonlength,1:latlength,1:altlength,j)
	999	else
	1000	out_4d_data(1:lonlength,1:latlength,1:altlength,j)=&
	1001	wk_4d_data(1:lonlength,2:latlength+1,1:altlength,j)
	1002	endif
	1003	enddo
	1004	! write interpolated data to output
	1005	corner(:)=1
	1006	edges(1)=lonlength
	1007	edges(2)=latlength
	1008	edges(3)=altlength
	1009	edges(4)=timelength
	1010	ierr=NF_PUT_VARA_REAL(outfid,var_id(i),corner,edges,out_4d_data)
	1011	if (ierr.ne.NF_NOERR) then
	1012	write(,) "Error: Could not write ",trim(var(i))," data to output file"
	1013	stop
	1014	else
	1015	write(,) "wrote ",trim(var(i))
	1016	endif
	1017	endif ! of if (tmpndims.eq.3)
	1018
	1019	enddo ! of do i=1,nbvar
	1020
	1021	!===============================================================================
	1022	! 6. Close output file
	1023	!===============================================================================
	1024	ierr=NF_CLOSE(outfid)
	1025	if (ierr.ne.NF_NOERR) then
	1026	write(,) 'Error, failed to close output file ',outfile
	1027	endif
	1028
	1029	end program hrecast
	1030
	1031
	1032
	1033	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1034	subroutine interp_horiz (varo,varn,imo,jmo,imn,jmn,lm, &
	1035	& rlonuo,rlatvo,rlonun,rlatvn)
	1036
	1037	!===========================================================
	1038	! Interpolation Horizontales des variables d'une grille LMDZ
	1039	! (des points SCALAIRES au point SCALAIRES)
	1040	! dans une autre grille LMDZ en conservant la quantite
	1041	! totale pour les variables intensives (/m2) : ex : Pression au sol
	1042	!
	1043	! Francois Forget (01/1995)
	1044	!===========================================================
	1045
	1046	IMPLICIT NONE
	1047
	1048	! Declarations:
	1049	! ==============
	1050	!
	1051	! ARGUMENTS
	1052	! """""""""
	1053
	1054	INTEGER imo, jmo ! dimensions ancienne grille (input)
	1055	INTEGER imn,jmn ! dimensions nouvelle grille (input)
	1056
	1057	REAL rlonuo(imo+1) ! Latitude et
	1058	REAL rlatvo(jmo) ! longitude des
	1059	REAL rlonun(imn+1) ! bord des
	1060	REAL rlatvn(jmn) ! cases "scalaires" (input)
	1061
	1062	INTEGER lm ! dimension verticale (input)
	1063	REAL varo (imo+1, jmo+1,lm) ! var dans l'ancienne grille (input)
	1064	REAL varn (imn+1,jmn+1,lm) ! var dans la nouvelle grille (output)
	1065
	1066	! Autres variables
	1067	! """"""""""""""""
	1068	INTEGER imnmx2,jmnmx2
	1069	! parameter (imnmx2=190,jmnmx2=100)
	1070	parameter (imnmx2=360,jmnmx2=190)
	1071	REAL airetest(imnmx2+1,jmnmx2+1)
	1072	INTEGER ii,jj,l
	1073
	1074	REAL airen ((imnmx2+1)*(jmnmx2+1)) ! aire dans la nouvelle grille
	1075	REAL airentotn ! aire totale pole nord dans la nouvelle grille
	1076	REAL airentots ! aire totale pole sud dans la nouvelle grille
	1077	! Info sur les ktotal intersection entre les cases new/old grille
	1078
	1079	! kmax: le nombre max d'intersections entre les 2 grilles horizontales
	1080	! On fixe kmax a la taille de la grille des donnees martiennes (360x179)
	1081	! + des pouiemes (cas ou une maille est a cheval sur 2 ou 4 mailles)
	1082	! Il y a un test dans iniinterp_h pour s'assurer que ktotal < kmax
	1083	INTEGER kmax, k, ktotal
	1084	parameter (kmax = 360*179 + 200000)
	1085	! parameter (kmax = 360*179 + 40000)
	1086
	1087	INTEGER iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
	1088	REAL intersec(kmax)
	1089	REAL R
	1090	REAL totn, tots
	1091	integer prev_sumdim
	1092	save prev_sumdim
	1093	data prev_sumdim /0/
	1094
	1095
	1096	logical firsttest, aire_ok
	1097	save firsttest
	1098	data firsttest /.true./
	1099	data aire_ok /.true./
	1100
	1101	integer imoS,jmoS,imnS,jmnS
	1102	save imoS,jmoS,imnS,jmnS
	1103	save ktotal,iik,jjk,jk,ik,intersec,airen
	1104	! REAL pi
	1105
	1106	! Test dimensions imnmx2 jmnmx2
	1107	!""""""""""""""""""""""""""""""
	1108	! test dimensionnement tableau airetest
	1109	if (imn.GT.imnmx2.OR.jmn.GT.jmnmx2) then
	1110	write(,) 'STOP pb dimensionnement tableau airetest'
	1111	write(,) 'il faut imn < imnmx2 et jmn < jmnmx2'
	1112	write(,) 'imn imnmx2', imn,imnmx2
	1113	write(,) 'jmn jmnmx2', jmn,jmnmx2
	1114	call exit(1)
	1115	endif
	1116
	1117	! initialisation
	1118	! --------------
	1119	! Si c'est le premier appel, on prepare l'interpolation
	1120	! en calculant pour chaque case autour d'un point scalaire de la
	1121	! nouvelle grille, la surface de intersection avec chaque
	1122	! case de l'ancienne grille.
	1123
	1124	! This must also be done if we change the dimension
	1125	if (imo+jmo+imn+jmn.ne.prev_sumdim) then
	1126	firsttest=.true.
	1127	prev_sumdim=imo+jmo+imn+jmn
	1128	end if
	1129
	1130	if (firsttest) then
	1131	call iniinterp_h(imo,jmo,imn,jmn ,kmax, &
	1132	& rlonuo,rlatvo,rlonun,rlatvn, &
	1133	& ktotal,iik,jjk,jk,ik,intersec,airen)
	1134	imoS=imo
	1135	jmoS=jmo
	1136	imnS=imn
	1137	jmnS=jmn
	1138	else
	1139	if(imo.NE.imoS.OR.jmo.NE.jmoS.OR.imn.NE.imnS.OR.jmn.NE.jmnS) then
	1140	call iniinterp_h(imo,jmo,imn,jmn ,kmax, &
	1141	& rlonuo,rlatvo,rlonun,rlatvn, &
	1142	& ktotal,iik,jjk,jk,ik,intersec,airen)
	1143	imoS=imo
	1144	jmoS=jmo
	1145	imnS=imn
	1146	jmnS=jmn
	1147	end if
	1148	end if
	1149
	1150	do l=1,lm
	1151	do jj =1 , jmn+1
	1152	do ii=1, imn+1
	1153	varn(ii,jj,l) =0.
	1154	end do
	1155	end do
	1156	end do
	1157
	1158	! Interpolation horizontale
	1159	! -------------------------
	1160	! boucle sur toute les ktotal intersections entre les cases
	1161	! de l'ancienne et la nouvelle grille
	1162	!
	1163
	1164	do k=1,ktotal
	1165	do l=1,lm
	1166	varn(iik(k),jjk(k),l) = varn(iik(k),jjk(k),l) &
	1167	& + varo(ik(k), jk(k),l)*intersec(k)/airen(iik(k) &
	1168	& +(jjk(k)-1)*(imn+1))
	1169	end do
	1170	end do
	1171
	1172	! Une seule valeur au pole pour les variables ! :
	1173	! -----------------------------------------------
	1174	DO l=1, lm
	1175	totn =0.
	1176	tots =0.
	1177
	1178
	1179	! moyenne du champ au poles (ponderee par les aires)
	1180	!"""""""""""""""""""""""""""""""
	1181	airentotn=0.
	1182	airentots=0.
	1183
	1184	do ii =1, imn+1
	1185	totn = totn + varn(ii,1,l)*airen(ii)
	1186	tots = tots + varn (ii,jmn+1,l)airen(jmn(imn+1)+ii)
	1187	airentotn=airentotn + airen(ii)
	1188	airentots=airentots + airen(jmn*(imn+1)+ii)
	1189	end do
	1190
	1191	do ii =1, imn+1
	1192	varn(ii,1,l) = totn/airentotn
	1193	varn(ii,jmn+1,l) = tots/airentots
	1194	end do
	1195
	1196	ENDDO
	1197
	1198
	1199	!---------------------------------------------------------------
	1200	! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
	1201	if (firsttest) then
	1202	! pi=2.*asin(1.)
	1203	firsttest = .false.
	1204	! write (,) 'INTERP. HORIZ. : TEST SUR LES AIRES:'
	1205
	1206	do jj =1 , jmn+1
	1207	do ii=1, imn+1
	1208	airetest(ii,jj) =0.
	1209	end do
	1210	end do
	1211	do k=1,ktotal
	1212	airetest(iik(k),jjk(k))= airetest(iik(k),jjk(k)) +intersec(k)
	1213	end do
	1214	do jj =1 , jmn+1
	1215	do ii=1, imn+1
	1216	r = airen(ii+(jj-1)*(imn+1))/airetest(ii,jj)
	1217	if ((r.gt.1.001).or.(r.lt.0.999)) then
	1218	write (,) '********** PROBLEME D'' AIRES !!!', &
	1219	& ' DANS L''INTERPOLATION HORIZONTALE'
	1220	write(,)'ii,jj,airen,airetest', &
	1221	& ii,jj,airen(ii+(jj-1)*(imn+1)),airetest(ii,jj)
	1222	aire_ok = .false.
	1223	end if
	1224	end do
	1225	end do
	1226	! if (aire_ok) write(,) 'INTERP. HORIZ. : AIRES OK'
	1227	endif
	1228
	1229	! FIN TEST FIN TEST FIN TEST FIN TEST FIN TEST FIN TEST FIN TEST
	1230	! --------------------------------------------------------------------
	1231
	1232
	1233	return
	1234	end
	1235	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1236	subroutine iniinterp_h (imo,jmo,imn,jmn ,kmax, &
	1237	& rlonuo,rlatvo,rlonun,rlatvn, &
	1238	& ktotal,iik,jjk,jk,ik,intersec,airen)
	1239
	1240	implicit none
	1241
	1242
	1243
	1244	! ---------------------------------------------------------
	1245	! Prepare l' interpolation des variables d'une grille LMDZ
	1246	! dans une autre grille LMDZ en conservant la quantite
	1247	! totale pour les variables intensives (/m2) : ex : Pression au sol
	1248	!
	1249	! (Pour chaque case autour d'un point scalaire de la nouvelle
	1250	! grille, on calcule la surface (en m2)en intersection avec chaque
	1251	! case de l'ancienne grille , pour la future interpolation)
	1252	!
	1253	! on calcule aussi l' aire dans la nouvelle grille
	1254	!
	1255	!
	1256	! Auteur: F.Forget 01/1995
	1257	! -------
	1258	!
	1259	! ---------------------------------------------------------
	1260	! Declarations:
	1261	! ==============
	1262	!
	1263	! ARGUMENTS
	1264	! """""""""
	1265	! INPUT
	1266	integer imo, jmo ! dimensions ancienne grille
	1267	integer imn,jmn ! dimensions nouvelle grille
	1268	integer kmax ! taille du tableau des intersections
	1269	real rlonuo(imo+1) ! Latitude et
	1270	real rlatvo(jmo) ! longitude des
	1271	real rlonun(imn+1) ! bord des
	1272	real rlatvn(jmn) ! cases "scalaires" (input)
	1273
	1274	! OUTPUT
	1275	integer ktotal ! nombre totale d''intersections reperees
	1276	integer iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
	1277	real intersec(kmax) ! surface des intersections (m2)
	1278	real airen (imn+1,jmn+1) ! aire dans la nouvelle grille
	1279
	1280
	1281
	1282
	1283	! Autres variables
	1284	! """"""""""""""""
	1285	integer i,j,ii,jj
	1286	integer imomx,jmomx,imnmx1,jmnmx1
	1287	! parameter (imomx=361,jmomx=180,imnmx1=190,jmnmx1=100)
	1288	parameter (imomx=361,jmomx=180,imnmx1=360,jmnmx1=190)
	1289	real a(imomx+1),b(imomx+1),c(jmomx+1),d(jmomx+1)
	1290	real an(imnmx1+1),bn(imnmx1+1)
	1291	real cn(jmnmx1+1),dn(jmnmx1+1)
	1292	real aa, bb,cc,dd
	1293	real pi
	1294
	1295	pi = 2.*ASIN( 1. )
	1296
	1297	! Test dimensions imnmx1 jmnmx1
	1298	!""""""""""""""""""""""""""""""
	1299	! test dimensionnement tableau airetest
	1300	if (imn.GT.imnmx1.OR.jmn.GT.jmnmx1) then
	1301	write(,) 'STOP pb dimensionnement'
	1302	write(,) 'il faut imn < imnmx1 et jmn < jmnmx1'
	1303	write(,) 'imn imnmx1', imn,imnmx1
	1304	write(,) 'jmn jmnmx1', jmn,jmnmx1
	1305	call exit(1)
	1306	endif
	1307
	1308	if (imo.GT.imomx.OR.jmo.GT.jmomx) then
	1309	write(,) 'STOP pb dimensionnement'
	1310	write(,) 'il faut imo < imomx et jmo < jmomx'
	1311	write(,) 'imo imomx', imo,imomx
	1312	write(,) 'jmo jmomx', jmo,jmomx
	1313	call exit(1)
	1314	endif
	1315
	1316	! On repere les frontieres des cases :
	1317	! ===================================
	1318	!
	1319	! Attention, on ruse avec des latitudes = 90 deg au pole.
	1320
	1321
	1322	! Ancienne grile
	1323	! """"""""""""""
	1324	a(1) = - rlonuo(imo+1)
	1325	b(1) = rlonuo(1)
	1326	do i=2,imo+1
	1327	a(i) = rlonuo(i-1)
	1328	b(i) = rlonuo(i)
	1329	end do
	1330
	1331	d(1) = pi/2
	1332	do j=1,jmo
	1333	c(j) = rlatvo(j)
	1334	d(j+1) = rlatvo(j)
	1335	end do
	1336	c(jmo+1) = -pi/2
	1337
	1338
	1339	! Nouvelle grille
	1340	! """""""""""""""
	1341	an(1) = - rlonun(imn+1)
	1342	bn(1) = rlonun(1)
	1343	do i=2,imn+1
	1344	an(i) = rlonun(i-1)
	1345	bn(i) = rlonun(i)
	1346	end do
	1347
	1348	dn(1) = pi/2
	1349	do j=1,jmn
	1350	cn(j) = rlatvn(j)
	1351	dn(j+1) = rlatvn(j)
	1352	end do
	1353	cn(jmn+1) = -pi/2
	1354
	1355	! Calcul de la surface des cases scalaires de la nouvelle grille
	1356	! ==============================================================
	1357	do ii=1,imn + 1
	1358	do jj = 1,jmn+1
	1359	airen(ii,jj) = (bn(ii)-an(ii))*(sin(dn(jj))-sin(cn(jj)))
	1360	end do
	1361	end do
	1362
	1363	! Calcul de la surface des intersections
	1364	! ======================================
	1365
	1366	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1367	! test dimenssion kmax (calcul de ktotal)
	1368	!"""""""""""""""""""""""""""""""""""""""
	1369	! Calcul de ktotal, mais ralonge beaucoup en temps => pour debug
	1370	! write(,)
	1371	! write(,) 'DEBUT DU TEST KMAX'
	1372	ktotal = 0
	1373	do jj = 1,jmn+1
	1374	do j=1, jmo+1
	1375	if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
	1376	do ii=1,imn + 1
	1377	do i=1, imo +1
	1378	if ( ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i))) &
	1379	& .or. ((an(ii).lt.b(i)-2pi).and.(bn(ii).gt.a(i)-2pi) &
	1380	& .and.(b(i)-2*pi.lt.-pi) ) &
	1381	& .or. ((an(ii).lt.b(i)+2pi).and.(bn(ii).gt.a(i)+2pi) &
	1382	& .and.(a(i)+2*pi.gt.pi) ) &
	1383	& )then
	1384	ktotal = ktotal +1
	1385	end if
	1386	enddo
	1387	enddo
	1388	end if
	1389	enddo
	1390	enddo
	1391
	1392	if (kmax.LT.ktotal) then
	1393	write(,)
	1394	write(,) '****** ATTENTION ******'
	1395	write(,) 'kmax =',kmax
	1396	write(,) 'ktotal =',ktotal
	1397	write(,) 'Changer la valeur de kmax dans interp_horiz.F '
	1398	write(,) 'avec kmax >= ktotal'
	1399	write(,) 'EXIT dans iniinterp_h'
	1400	call exit(1)
	1401	else
	1402	! write(,) 'kmax =',kmax
	1403	! write(,) 'ktotal =',ktotal
	1404	end if
	1405	! write(,) 'FIN DU TEST KMAX'
	1406	! write(,)
	1407	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1408
	1409	! boucle sur la nouvelle grille
	1410	! """"""""""""""""""""""""""""
	1411	ktotal = 0
	1412	do jj = 1,jmn+1
	1413	do j=1, jmo+1
	1414	if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
	1415	do ii=1,imn + 1
	1416	do i=1, imo +1
	1417	if ( ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i))) &
	1418	& .or. ((an(ii).lt.b(i)-2pi).and.(bn(ii).gt.a(i)-2pi) &
	1419	& .and.(b(i)-2*pi.lt.-pi) ) &
	1420	& .or. ((an(ii).lt.b(i)+2pi).and.(bn(ii).gt.a(i)+2pi) &
	1421	& .and.(a(i)+2*pi.gt.pi) ) &
	1422	& )then
	1423	ktotal = ktotal +1
	1424	iik(ktotal) =ii
	1425	jjk(ktotal) =jj
	1426	ik(ktotal) =i
	1427	jk(ktotal) =j
	1428
	1429	dd = min(d(j), dn(jj))
	1430	cc = cn(jj)
	1431	if (cc.lt. c(j))cc=c(j)
	1432	if((an(ii).lt.b(i)-2*pi).and. &
	1433	& (bn(ii).gt.a(i)-2*pi)) then
	1434	bb = min(b(i)-2*pi,bn(ii))
	1435	aa = an(ii)
	1436	if (aa.lt.a(i)-2pi) aa=a(i)-2pi
	1437	else if((an(ii).lt.b(i)+2*pi).and. &
	1438	& (bn(ii).gt.a(i)+2*pi)) then
	1439	bb = min(b(i)+2*pi,bn(ii))
	1440	aa = an(ii)
	1441	if (aa.lt.a(i)+2pi) aa=a(i)+2pi
	1442	else
	1443	bb = min(b(i),bn(ii))
	1444	aa = an(ii)
	1445	if (aa.lt.a(i)) aa=a(i)
	1446	end if
	1447	intersec(ktotal)=(bb-aa)*(sin(dd)-sin(cc))
	1448	end if
	1449	end do
	1450	end do
	1451	end if
	1452	end do
	1453	end do
	1454
	1455
	1456	! TEST INFO
	1457	! DO k=1,ktotal
	1458	! ii = iik(k)
	1459	! jj = jjk(k)
	1460	! i = ik(k)
	1461	! j = jk(k)
	1462	! if ((ii.eq.10).and.(jj.eq.10).and.(i.eq.10).and.(j.eq.10))then
	1463	! if (jj.eq.2.and.(ii.eq.1))then
	1464	! write(,) '**************** jj=',jj,'ii=',ii
	1465	! write(,) 'i,j =',i,j
	1466	! write(,) 'an,bn,cn,dn', an(ii), bn(ii), cn(jj),dn(jj)
	1467	! write(,) 'a,b,c,d', a(i), b(i), c(j),d(j)
	1468	! write(,) 'intersec(k)',intersec(k)
	1469	! write(,) 'airen(ii,jj)=',airen(ii,jj)
	1470	! end if
	1471	! END DO
	1472
	1473	return
	1474	end

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

Download in other formats:

Original Format