Context Navigation

api.F90 @ 73

Last change on this file since 73 was 73, checked in by aslmd, 14 years ago
LMD_MM_MARS: corrections cycle de l'eau propagees a la nouvelle physique. + corrections readmeteo.F90 [version synchronisee precedemment n etait pas la plus a jour] + corrections api.F90 pour avoir cp, R comme GCM
File size: 69.4 KB

Rev	Line
[19]	1	PROGRAM api
	2
	3	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	4	! API !
	5	! Altitude and Pressure Interpolator !
	6	! !
	7	! This program is based on !
	8	! ------------------------ !
	9	! p_interp v1.0 !
	10	! http://www.mmm.ucar.edu/wrf/src/p_interp.tar.gz !
	11	! September 2008 - Cindy Bruyere [NCAR, USA] !
	12	! !
	13	! Modifications !
	14	! ------------- !
	15	! - Presentation of the program (cosmetics) !
	16	! - Additional routine and arguments in existing routines !
	17	! - Improve memory management !
	18	! - Martian constants + Addition of 'grav' variable !
	19	! - Interpolation to height above areoid (MOLA zero datum) !
	20	! October 2009 - Aymeric Spiga [The Open University, UK] !
	21	! - Interpolation to height above surface !
	22	! - Rotated winds (e.g. for polar projections) !
	23	! November 2009 - AS !
	24	! !
	25	! Purpose !
	26	! ------- !
	27	! Program to read wrfout data (possibly several files) !
	28	! and interpolate to pressure or height levels !
	29	! A new NETCDF file is generated with appropriate suffix !
	30	! The program reads namelist.api !
	31	! !
	32	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	33
	34	IMPLICIT NONE
	35	INCLUDE 'netcdf.inc'
	36
	37	!!
	38	!! EARTH CONSTANTS
	39	!!
	40	!REAL, PARAMETER :: Rd = 287.04 ! gas constant m2 s-2 K-1
	41	!REAL, PARAMETER :: Cp = 7.Rd/2. ! r=8.314511E+01000.E+0/mugaz
	42	!REAL, PARAMETER :: RCP = Rd/Cp
	43	!REAL, PARAMETER :: p0 = 100000.
	44	!REAL, PARAMETER :: grav = 9.81
	45	!REAL, PARAMETER :: tpot0 = 300.
	46
	47	!
	48	! MARS CONSTANTS
	49	!
[73]	50	!REAL, PARAMETER :: Rd = 192. ! gas constant m2 s-2 K-1
	51	!REAL, PARAMETER :: Cp = 844.6 ! r= 8.314511E+0*1000.E+0/mugaz
	52	REAL, PARAMETER :: Rd = 191.0
	53	REAL, PARAMETER :: Cp = 744.5
[19]	54	REAL, PARAMETER :: RCP = Rd/Cp
	55	REAL, PARAMETER :: p0 = 610.
	56	REAL, PARAMETER :: grav = 3.72
	57	REAL, PARAMETER :: tpot0 = 220. ! reference potential temperature
	58
	59	!
	60	! VARIABLES
	61	!
	62	CHARACTER, ALLOCATABLE, DIMENSION(:,:,:,:) :: text
	63	CHARACTER (LEN=31),ALLOCATABLE, DIMENSION(:) :: dnamei, dnamej
	64	CHARACTER(LEN=250),ALLOCATABLE, DIMENSION(:) :: input_file_names
	65	CHARACTER(LEN=250),ALLOCATABLE, DIMENSION(:) :: output_file_names
	66	DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:) :: ddata1, ddata2
	67	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: data1, data2, data3
	68	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: pres_field, pres_out
	69	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: pres_stagU, pres_stagV
	70	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: ght, phb, qv, tk, rh, tpot, u, v, umet, vmet
	71	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: psfc
	72	REAL, ALLOCATABLE, DIMENSION(:,:) :: ter
	73	REAL, ALLOCATABLE, DIMENSION(:,:) :: longi, lati
	74	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: longit, latit
	75	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: interm1, interm2
	76	INTEGER, ALLOCATABLE, DIMENSION(:) :: dvali, dvalj
	77	INTEGER, ALLOCATABLE, DIMENSION(:,:,:,:) :: idata1, idata2
	78	INTEGER, DIMENSION(4) :: start_dims = 1
	79	INTEGER, DIMENSION(4) :: dims_in, dims_out
	80	INTEGER, DIMENSION(6) :: ishape, jshape
	81	CHARACTER (LEN=80) :: cval
	82	CHARACTER (LEN=31) :: cname, test_dim_name
	83	CHARACTER (LEN=80) :: input_file, output_file, att_text
	84	CHARACTER (LEN=250) :: path_to_input
	85	CHARACTER (LEN=250) :: path_to_output
	86	CHARACTER (LEN=250) :: input_name
	87	CHARACTER (LEN=250) :: output_name, tmp_name
	88	CHARACTER (LEN=10) :: option
	89	CHARACTER (LEN=132) :: command
	90	CHARACTER (LEN=20) :: process, dummy
	91	CHARACTER (LEN=2000) :: fields, process_these_fields
	92	REAL, DIMENSION(99) :: interp_levels
	93	REAL :: rval
	94	REAL :: MISSING=1.e36
	95	REAL :: truelat1, truelat2, stand_lon
	96	INTEGER :: map_proj
	97	INTEGER :: LINLOG = 1
	98	INTEGER :: interp_method=1
	99	INTEGER :: extrapolate=0
	100	INTEGER :: ncid, mcid, rcode
	101	INTEGER :: idm, ndims, nvars, natt, ngatts
	102	INTEGER :: nunlimdimid
	103	INTEGER :: i, ii, j, jj, ix, iy, iweg, isng, ibtg
	104	INTEGER :: ivar, jvar
	105	INTEGER :: times_in_file
	106	INTEGER :: ilen, itype, ival, na, funit, ios
	107	INTEGER :: num_metgrid_levels
	108	INTEGER :: number_of_input_files
	109	INTEGER :: ierr, loop, loslen, strlen, lent
	110	INTEGER :: is_there
	111	INTEGER :: kk
	112	LOGICAL :: is_used
	113	LOGICAL :: debug=.FALSE.
	114	LOGICAL :: interpolate=.FALSE.
	115	LOGICAL :: unstagger_grid=.FALSE.
	116	LOGICAL :: fix_meta_stag=.FALSE.
	117	LOGICAL :: bit64=.FALSE.
	118	LOGICAL :: first=.TRUE.
	119	LOGICAL :: oldvar=.FALSE.
	120
	121	!
	122	! NAMELISTS
	123	!
	124	NAMELIST /io/ path_to_input, input_name, path_to_output, output_name, &
	125	process, fields, debug, bit64, oldvar
	126	NAMELIST /interp_in/ interp_levels, interp_method, extrapolate, unstagger_grid
	127
	128	!
	129	! DEFAULT VALUES for VARIABLES
	130	!
	131	path_to_input = './'
	132	path_to_output = './'
	133	output_name = ' '
	134	interp_levels = -99999.
	135	process = 'all'
	136
	137	!
	138	! READ NAMELIST
	139	!
	140	DO funit=10,100
	141	INQUIRE(unit=funit, opened=is_used)
	142	IF (.not. is_used) EXIT
	143	END DO
	144	OPEN(funit,file='namelist.api',status='old',form='formatted',iostat=ios)
	145	IF ( ios /= 0 ) STOP "ERROR opening namelist.api"
	146	READ(funit,io)
	147	READ(funit,interp_in)
	148	CLOSE(funit)
	149
	150	!
	151	! INPUT FILE NAMES
	152	!
	153	lent = len_trim(path_to_input)
	154	IF ( path_to_input(lent:lent) /= "/" ) THEN
	155	path_to_input = TRIM(path_to_input)//"/"
	156	ENDIF
	157	lent = len_trim(path_to_output)
	158	IF ( path_to_output(lent:lent) /= "/" ) THEN
	159	path_to_output = TRIM(path_to_output)//"/"
	160	ENDIF
	161	input_name = TRIM(path_to_input)//TRIM(input_name)
	162	!
	163	! BUILD a UNIX COMMAND based on "ls" of all of the input files
	164	!
	165	loslen = LEN ( command )
	166	CALL all_spaces ( command , loslen )
	167	WRITE ( command , FMT='("ls -1 ",A," > .foo")' ) TRIM ( input_name )
	168
	169	!
	170	! NUMBER OF INPUTS FILES
	171	!
	172	! We stuck all of the matching files in the ".foo" file. Now we place the
	173	! number of the those file (i.e. how many there are) in ".foo1".
	174	CALL SYSTEM ( TRIM ( command ) )
	175	CALL SYSTEM ( '( cat .foo \| wc -l > .foo1 )' )
	176	! Read the number of files.
	177	OPEN (FILE = '.foo1' , &
	178	UNIT = 112 , &
	179	STATUS = 'OLD' , &
	180	ACCESS = 'SEQUENTIAL' , &
	181	FORM = 'FORMATTED' )
	182	READ ( 112 , * ) number_of_input_files
	183	CLOSE ( 112 )
	184	! If there are zero files, we are toast.
	185	IF ( number_of_input_files .LE. 0 ) THEN
	186	print*, ' Oops, we need at least ONE input file for the program to read.'
	187	print*, ' Make sure you have the path, and name file(s) correct,'
	188	print*, ' including wild characters if needed.'
	189	STOP
	190	END IF
	191
	192	!
	193	! GET FILENAMES IN THE PROGRAM
	194	!
	195	! Allocate space for this many files.
	196	ALLOCATE ( input_file_names(number_of_input_files) , STAT=ierr )
	197	ALLOCATE ( output_file_names(number_of_input_files) , STAT=ierr )
	198	! Did the allocate work OK?
	199	IF ( ierr .NE. 0 ) THEN
	200	print*, ' tried to allocate ', number_of_input_files, ' input files, (look at ./foo)'
	201	STOP
	202	END IF
	203	! Initialize all of the file names to blank.
	204	input_file_names = ' ' // &
	205	' ' // &
	206	' '
	207	output_file_names = ' ' // &
	208	' ' // &
	209	' '
	210	! Open the file that has the list of filenames.
	211	OPEN (FILE = '.foo' , &
	212	UNIT = 111 , &
	213	STATUS = 'OLD' , &
	214	ACCESS = 'SEQUENTIAL' , &
	215	FORM = 'FORMATTED' )
	216	! Read all of the file names and store them - define output file name.
	217	LINLOG = interp_method
	218	DO loop = 1 , number_of_input_files
	219	READ ( 111 , FMT='(A)' ) input_file_names(loop)
	220	IF ( output_name == ' ' ) THEN
	221	ilen = INDEX(TRIM(input_file_names(loop)),'/',.TRUE.)
	222	output_file_names(loop) = TRIM(path_to_output)//input_file_names(loop)(ilen+1:)
	223	IF ( LINLOG .le. 2 ) output_file_names(loop) = TRIM(output_file_names(loop))//"_p"
	224	IF ( LINLOG == 3 ) output_file_names(loop) = TRIM(output_file_names(loop))//"_z"
	225	IF ( LINLOG == 4 ) output_file_names(loop) = TRIM(output_file_names(loop))//"_zabg"
	226	ELSE
	227	IF ( number_of_input_files == 1 ) THEN
	228	output_file_names(loop) = TRIM(path_to_output)//TRIM(output_name)
	229	ELSE
	230	write(tmp_name,'(A,A,"_",I4.4)') TRIM(path_to_output), TRIM(output_name), loop
	231	output_file_names(loop) = tmp_name
	232	ENDIF
	233	ENDIF
	234	END DO
	235	CLOSE ( 112 )
	236	print*, " "
	237	! We clean up our own messes.
	238	CALL SYSTEM ( '/bin/rm -f .foo' )
	239	CALL SYSTEM ( '/bin/rm -f .foo1' )
	240
	241	!
	242	! LIST OF FIELD WANTED for OUTPUT
	243	!
	244	! Do we have a list of field that we want on output?
	245	! process_these_fields = ','
	246	!!!! ICI CHAMPS SORTIS PAR DEFAULT
	247	!!!! ICI CHAMPS SORTIS PAR DEFAULT
	248	!!!! ICI CHAMPS SORTIS PAR DEFAULT
	249	process_these_fields = ',Times,XLAT,XLONG,' !! les premiere et derniere virgules sont importantes
	250	!!!! ICI CHAMPS SORTIS PAR DEFAULT
	251	!!!! ICI CHAMPS SORTIS PAR DEFAULT
	252	!!!! ICI CHAMPS SORTIS PAR DEFAULT
	253	IF ( INDEX(process,'list') /= 0) THEN
	254	DO i = 1 , len(fields)
	255	IF (fields(i:i) /= ' ' ) THEN
	256	process_these_fields = trim(process_these_fields)//fields(i:i)
	257	ENDIF
	258	END DO
	259	process_these_fields = trim(process_these_fields)//","
	260	END IF
	261
	262	!
	263	! HELLO to WORLD
	264	!
	265	write(6,*)
	266	write(6,*) "##############################################################"
	267	write(6,'(A,i4,A)') " RUNNING api on ", number_of_input_files, " file(s)."
	268	write(6,*)
	269	write(6,'(A,$)') " A) INTERPOLATION METHOD: "
	270	IF ( LINLOG == 1 ) write(6,*) " Pressure - linear in p"
	271	IF ( LINLOG == 2 ) write(6,*) " Pressure - linear in log p"
	272	IF ( LINLOG == 3 ) write(6,*) " Height above areoid (MOLA zero datum)"
	273	IF ( LINLOG == 4 ) write(6,*) " Height above surface"
	274	write(6,'(A,$)') " B) VERTICAL EXTRAPOLATION: "
	275	IF (extrapolate == 0) write(6,*) " BELOW GROUND and ABOVE model top will be set to missing values "
	276	IF (extrapolate == 1) write(6,*) " BELOW GROUND will be extrapolated and ABOVE model top will be set to values at model top"
	277	write(6,'(A,$)') " C) HORIZONTAL GRID: "
	278	IF (.not. unstagger_grid) write(6,*) " Data will be output on C-grid "
	279	IF (unstagger_grid) write(6,*) " Data will be output on unstaggered grid "
	280
	281	!
	282	! GET LEVELS (pressure or altitude) TO INTERPOLATE TO
	283	!
	284	write(6,*)
	285	IF ( LINLOG .le. 2 ) write(6,*) "INTERPOLATING TO PRESSURE LEVELS: "
	286	IF ( LINLOG == 3 ) write(6,*) "INTERPOLATING TO ALTITUDE LEVELS: "
	287	IF ( LINLOG == 4 ) write(6,*) "INTERPOLATING TO ABG LEVELS: "
	288	num_metgrid_levels = 0
	289	DO
	290	IF (interp_levels(num_metgrid_levels+1) == -99999.) EXIT
	291	num_metgrid_levels = num_metgrid_levels + 1
	292	if (mod(num_metgrid_levels,8) /= 0 )write(6,'(f8.3,$)') interp_levels(num_metgrid_levels)
	293	if (mod(num_metgrid_levels,8) == 0 )write(6,'(f8.3)') interp_levels(num_metgrid_levels)
	294	IF ( LINLOG .le. 2 ) interp_levels(num_metgrid_levels) = interp_levels(num_metgrid_levels) * 100.0 !!! Pa
	295	IF ( LINLOG .gt. 2 ) interp_levels(num_metgrid_levels) = interp_levels(num_metgrid_levels) * 1000.0 !!! km
	296	END DO
	297	write(6,*)
	298	write(6,*)
	299
	300	!
	301	! LOOP on FILES
	302	!
	303	DO loop = 1, number_of_input_files
	304
	305	IF (debug) write(6,*) "##############################################################"
	306	input_file = input_file_names(loop)
	307	output_file = output_file_names(loop)
	308
	309	IF ( .not. debug ) write(6,*) " Output will be written to: ",trim(output_file)
	310
	311	IF (debug) THEN
	312	write(6,*) " INPUT FILE: ",trim(input_file)
	313	write(6,*) " OUTPUT FILE: ",trim(output_file)
	314	write(6,*) " "
	315	ENDIF
	316
	317	!
	318	! OPEN INPUT AND OUTPUT FILE
	319	!
	320	rcode = nf_open(input_file, 0, ncid)
	321	if (rcode .ne. nf_noerr) call handle_err(rcode)
	322	if (bit64) then
	323	rcode = nf_create(output_file, NF_64BIT_OFFSET, mcid)
	324	else
	325	rcode = nf_create(output_file, 0, mcid)
	326	endif
	327	if (rcode .ne. nf_noerr) call handle_err(rcode)
	328
	329	!
	330	! GET BASIC INFORMATION ABOUT THE FILE
	331	! most important
	332	! ndims: number of dimensions
	333	! nvars: number of variables
	334	! ngatts: number of global attributes
	335	rcode = nf_inq(ncid, ndims, nvars, ngatts, nunlimdimid)
	336	if (rcode .ne. nf_noerr) call handle_err(rcode)
	337	IF (debug) THEN
	338	write(6,*) ' INPUT file has = ',ndims, ' dimensions, '
	339	write(6,*) ' ',nvars, ' variables, and '
	340	write(6,*) ' ',ngatts,' global attributes '
	341	write(6,*) " "
	342	ENDIF
	343	rcode = nf_get_att_int (ncid, nf_global, 'WEST-EAST_GRID_DIMENSION', iweg)
	344	rcode = nf_get_att_int (ncid, nf_global, 'SOUTH-NORTH_GRID_DIMENSION', isng)
	345	rcode = nf_get_att_int (ncid, nf_global, 'BOTTOM-TOP_GRID_DIMENSION', ibtg)
	346
	347	!
	348	! ALLOCATE SOME VARIABLES
	349	!
	350	IF (ALLOCATED(dnamei)) deallocate(dnamei)
	351	ALLOCATE (dnamei(20))
	352	IF (ALLOCATED(dnamej)) deallocate(dnamej)
	353	ALLOCATE (dnamej(20))
	354	IF (ALLOCATED(dvali)) deallocate(dvali)
	355	ALLOCATE (dvali(20))
	356	IF (ALLOCATED(dvalj)) deallocate(dvalj)
	357	ALLOCATE (dvalj(20))
	358
	359	!
	360	! READ ALL DIMS FROM INPUT FILE AND CREATE SOME DIMS FOR OUTPUT FILE
	361	!
	362	j = 0
	363	DO i = 1, ndims
	364	rcode = nf_inq_dim(ncid, i, dnamei(i), dvali(i))
	365	IF ( dnamei(i) == "Time" ) THEN
	366	j = j + 1
	367	dnamej(j) = dnamei(i)
	368	dvalj(j) = dvali(i)
	369	rcode = nf_def_dim(mcid, dnamej(j), NF_UNLIMITED, j)
	370	times_in_file = dvali(i)
	371	ENDIF
	372	ENDDO
	373	!!! Create the new height/pressure dims
	374	j = j + 1
	375	IF ( LINLOG .le. 2 ) THEN
	376	dnamej(j) = 'pressure'
	377	dvalj(j) = num_metgrid_levels
	378	ELSE IF ( LINLOG .eq. 3 ) THEN
	379	dnamej(j) = 'altitude'
	380	dvalj(j) = num_metgrid_levels
	381	ELSE
	382	dnamej(j) = 'altitude_abg'
	383	dvalj(j) = num_metgrid_levels
	384	ENDIF
	385	rcode = nf_def_dim(mcid, dnamej(j), dvalj(j), j)
	386
	387	!
	388	! DEALING WITH THE GLOBAL ATTRIBUTES
	389	!
	390	IF (debug) THEN
	391	write(6,*)
	392	write(6,*) " OUTPUT FILE attributes:"
	393	ENDIF
	394	do i = 1, ngatts
	395	rcode = nf_inq_attname(ncid, nf_global, i, cname)
	396	rcode = nf_inq_atttype(ncid, nf_global, cname, itype)
	397	rcode = nf_inq_attlen (ncid, nf_global, cname, ilen)
	398	if ( itype .eq. 2 ) then ! characters
	399	rcode = nf_get_att_text (ncid, nf_global, cname, cval)
	400	if(cname(1:5) .eq. 'TITLE') then
	401	IF ( LINLOG .le. 2 ) THEN
	402	cval = cval(1:ilen)//" - ON PRES LEVELS"
	403	ELSE
	404	cval = cval(1:ilen)//" - ON z LEVELS"
	405	ENDIF
	406	ilen = len_trim(cval)
	407	endif
	408	IF (debug) &
	409	write(6,'(" i = ",i2," : ",A," = ",A)') &
	410	i,cname,cval(1:ilen)
	411	rcode = nf_put_att_text(mcid, nf_global, cname, ilen,&
	412	cval(1:ilen))
	413	elseif ( itype .eq. 4 ) then ! integers
	414	rcode = nf_get_att_int (ncid, nf_global, cname, ival)
	415	IF ( INDEX(cname,'BOTTOM-TOP_PATCH') == 0 ) THEN
	416	IF (cname .eq. 'BOTTOM-TOP_GRID_DIMENSION') ival = num_metgrid_levels
	417	IF (debug) &
	418	write(6,'(" i = ",i2," : ",A," = ",i7)') &
	419	i,cname,ival
	420	rcode = nf_put_att_int(mcid, nf_global, cname, itype,&
	421	ilen, ival)
	422	ENDIF
	423	IF (cname .eq. 'MAP_PROJ') map_proj = ival
	424	elseif ( itype .eq. 5 ) then ! real
	425	rcode = nf_get_att_real (ncid, nf_global, cname, rval)
	426	IF (debug) &
	427	write(6,'(" i = ",i2," : ",A," = ",G18.10E2)') &
	428	i,cname,rval
	429	rcode = nf_put_att_real(mcid, nf_global, cname, itype,&
	430	ilen, rval)
	431	IF (cname .eq. 'TRUELAT1') truelat1 = rval
	432	IF (cname .eq. 'TRUELAT2') truelat2 = rval
	433	IF (cname .eq. 'STAND_LON') stand_lon = rval
	434	end if
	435	enddo
	436	rcode = nf_enddef(mcid)
	437
	438	!
	439	! WE NEED SOME BASIC FIELDS
	440	!
	441	! --> P_TOP
	442	!
	443	IF ( LINLOG .le. 2 ) THEN
	444	IF ( DEBUG ) PRINT *, 'P_TOP'
	445	IF (ALLOCATED(data1)) deallocate(data1)
	446	allocate (data1(times_in_file,1,1,1))
	447	rcode = nf_inq_varid ( ncid, "P_TOP", i )
	448	rcode = nf_get_var_real ( ncid, i, data1 )
	449	IF ( first ) THEN
	450	IF ( extrapolate == 1 .AND. &
	451	(data1(1,1,1,1)-interp_levels(num_metgrid_levels)) > 0.0 ) THEN
	452	write(6,*)
	453	write(6,*) "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
	454	write(6,*) " WARNING: Highest requested pressure level is above PTOP."
	455	write(6,'(A,F7.2,A)') " Use all pressure level data above", data1(1,1,1,1)*.01, " mb"
	456	write(6,*) " with caution."
	457	write(6,*) "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
	458	ENDIF
	459	first = .FALSE.
	460	ENDIF
	461	deallocate (data1)
	462	ENDIF
	463	!
	464	! --> INTERPOLATION LEVELS
	465	!
	466	IF (ALLOCATED(pres_out)) deallocate(pres_out)
	467	allocate (pres_out(iweg-1, isng-1, num_metgrid_levels, times_in_file))
	468	do i = 1, num_metgrid_levels
	469	pres_out (:,:,i,:) = interp_levels(i)
	470	enddo
	471	!
	472	! --> LONGITUDES + everything for ROTATED WINDS
	473	!
	474	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'uvmet') /= 0 ) THEN
	475	IF (ALLOCATED(longi)) deallocate(longi)
	476	IF (ALLOCATED(longit)) deallocate(longit)
	477	allocate (longi(iweg-1, isng-1))
	478	allocate (longit(iweg-1, isng-1, times_in_file))
	479	IF ( DEBUG ) PRINT *, 'LONGITUDE'
	480	rcode = nf_inq_varid ( ncid, "XLONG", i )
	481	rcode = nf_get_var_real ( ncid, i, longit )
	482	longi = longit(:,:,1)
	483	deallocate(longit)
	484	IF (ALLOCATED(lati)) deallocate(lati)
	485	IF (ALLOCATED(latit)) deallocate(latit)
	486	allocate (lati(iweg-1, isng-1))
	487	allocate (latit(iweg-1, isng-1, times_in_file))
	488	IF ( DEBUG ) PRINT *, 'LATITUDE'
	489	rcode = nf_inq_varid ( ncid, "XLAT", i )
	490	rcode = nf_get_var_real ( ncid, i, latit )
	491	lati = latit(:,:,1)
	492	deallocate(latit)
	493	IF (ALLOCATED(u)) deallocate(u)
	494	IF (oldvar) THEN
	495	allocate (u (iweg, isng-1, ibtg-1, times_in_file ))
	496	IF ( DEBUG ) PRINT *, 'U COMPONENT'
	497	rcode = nf_inq_varid ( ncid, "U", i )
	498	rcode = nf_get_var_real ( ncid, i, u )
	499	ELSE
	500	allocate (u (iweg-1, isng-1, ibtg-1, times_in_file ))
	501	IF ( DEBUG ) PRINT *, 'UAVE COMPONENT'
	502	rcode = nf_inq_varid ( ncid, "UAVE", i )
	503	rcode = nf_get_var_real ( ncid, i, u )
	504	ENDIF
	505	IF (ALLOCATED(v)) deallocate(v)
	506	IF (oldvar) THEN
	507	allocate (v (iweg-1, isng, ibtg-1, times_in_file ))
	508	IF ( DEBUG ) PRINT *, 'V COMPONENT'
	509	rcode = nf_inq_varid ( ncid, "V", i )
	510	rcode = nf_get_var_real ( ncid, i, v )
	511	ELSE
	512	allocate (v (iweg-1, isng-1, ibtg-1, times_in_file ))
	513	IF ( DEBUG ) PRINT *, 'VAVE COMPONENT'
	514	rcode = nf_inq_varid ( ncid, "VAVE", i )
	515	rcode = nf_get_var_real ( ncid, i, v )
	516	ENDIF
	517	IF (ALLOCATED(umet)) deallocate(umet)
	518	allocate (umet (iweg-1, isng-1, ibtg-1, times_in_file ))
	519	IF (ALLOCATED(vmet)) deallocate(vmet)
	520	allocate (vmet (iweg-1, isng-1, ibtg-1, times_in_file ))
	521	IF (ALLOCATED(interm1)) deallocate(interm1)
	522	allocate (interm1(iweg-1, isng-1, ibtg-1 ))
	523	IF (ALLOCATED(interm2)) deallocate(interm2)
	524	allocate (interm2(iweg-1, isng-1, ibtg-1 ))
	525	IF ( DEBUG ) PRINT *, 'ROTATE WINDS'
	526	write(6,*) " Data will be output on unstaggered grid "
	527	do kk = 1, times_in_file
	528	IF ( DEBUG ) print *, kk
	529	IF (oldvar) THEN
	530	interm1(1:iweg-1,:,:) = ( u(1:iweg-1,:,:,kk) + u(2:iweg,:,:,kk) ) * .5
	531	interm2(:,1:isng-1,:) = ( v(:,1:isng-1,:,kk) + v(:,2:isng,:,kk) ) * .5
	532	ELSE
	533	interm1(1:iweg-1,:,:) = u(1:iweg-1,:,:,kk)
	534	interm2(:,1:isng-1,:) = v(:,1:isng-1,:,kk)
	535	ENDIF
	536	if (unstagger_grid) map_proj=3 !!! AS: unstagger_grid=T makes the program to put unstaggered U and V in uvmet
	537	CALL calc_uvmet( interm1, interm2, umet(:,:,:,kk), vmet(:,:,:,kk), &
	538	truelat1, truelat2, stand_lon, map_proj, longi, lati, iweg-1, isng-1, ibtg-1 )
	539	enddo
	540	deallocate(lati)
	541	deallocate(longi)
	542	deallocate(u)
	543	deallocate(v)
	544	deallocate(interm1)
	545	deallocate(interm2)
	546	ENDIF
	547
	548	!
	549	! --> GEOPOTENTIAL HEIGHT (in mode 1-2, no need if geopotential is not requested)
	550	!
	551	IF ( LINLOG .gt. 2 .OR. INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'GHT') /= 0 ) THEN
	552	IF ( DEBUG ) PRINT *, 'GEOPOTENTIAL'
	553	!IF (ALLOCATED(ght)) deallocate(ght)
	554	!allocate (ght(iweg-1, isng-1, ibtg-1, times_in_file))
	555	IF (ALLOCATED(phb)) deallocate(phb)
	556	allocate (phb(iweg-1, isng-1, ibtg, times_in_file ))
	557	! IF (ALLOCATED(data1)) deallocate(data1)
	558	! allocate (data1(iweg-1, isng-1, ibtg, times_in_file ))
	559	! rcode = nf_inq_varid ( ncid, "PH", i )
	560	! rcode = nf_get_var_real ( ncid, i, data1 )
	561	! rcode = nf_inq_varid ( ncid, "PHB", i )
	562	! rcode = nf_get_var_real ( ncid, i, phb )
	563	! data1 = (data1 + phb)
	564	! ght(:,:,1:ibtg-1,:) = ( data1(:,:,1:ibtg-1,:) + data1(:,:,2:ibtg,:) )*.5
	565	! deallocate (data1)
	566	rcode = nf_inq_varid ( ncid, "PHTOT", i )
	567	rcode = nf_get_var_real ( ncid, i, phb )
	568	!ght(:,:,1:ibtg-1,:) = phb(:,:,1:ibtg-1,:) !! costly in memory !!
	569	!deallocate (phb)
	570	ENDIF
	571	!
	572	! --> PRESSURE (in mode 3-4, no need if regular temperature is not requested)
	573	!
	574	IF ( LINLOG .le. 2 .OR. INDEX(process,'all') /= 0 &
	575	.OR. INDEX(process_these_fields,'tk') /= 0 &
	576	.OR. INDEX(process_these_fields,'tpot') /= 0 ) THEN
	577	IF ( DEBUG ) PRINT *, 'PRESSURE'
	578	IF (ALLOCATED(pres_field)) deallocate(pres_field)
	579	allocate (pres_field(iweg-1, isng-1, ibtg-1, times_in_file ))
	580	! IF (ALLOCATED(data1)) deallocate(data1)
	581	! allocate (data1(iweg-1, isng-1, ibtg-1, times_in_file ))
	582	! rcode = nf_inq_varid ( ncid, "P", i )
	583	! rcode = nf_get_var_real ( ncid, i, pres_field )
	584	! rcode = nf_inq_varid ( ncid, "PB", i )
	585	! rcode = nf_get_var_real ( ncid, i, data1 )
	586	! pres_field = pres_field + data1
	587	rcode = nf_inq_varid ( ncid, "PTOT", i )
	588	rcode = nf_get_var_real ( ncid, i, pres_field )
	589	! deallocate (data1)
	590	!
	591	! --> REGULAR and POTENTIAL TEMPERATURE
	592	!
	593	IF ( DEBUG ) PRINT *, 'TEMP'
	594	IF (ALLOCATED(tpot)) deallocate(tpot)
	595	allocate (tpot(iweg-1, isng-1, ibtg-1, times_in_file ))
	596	IF (oldvar) THEN
	597	rcode = nf_inq_varid ( ncid, "T", i )
	598	rcode = nf_get_var_real ( ncid, i, tpot )
	599	tpot = tpot + tpot0
	600	ELSE
	601	rcode = nf_inq_varid ( ncid, "TAVE", i )
	602	rcode = nf_get_var_real ( ncid, i, tpot )
	603	ENDIF
	604	IF ( LINLOG .le. 2 .OR. INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'tk') /= 0 ) THEN
	605	IF (ALLOCATED(tk)) deallocate(tk)
	606	allocate (tk(iweg-1, isng-1, ibtg-1, times_in_file ))
	607	tk = tpot * ( pres_field / p0 )**RCP
	608	ENDIF
	609	IF (INDEX(process_these_fields,'tpot') == 0 .AND. INDEX(process,'all') == 0) deallocate(tpot)
	610	ENDIF
	611	!
	612	! --> SURFACE PRESSURE
	613	!
	614	IF ( LINLOG .le. 2 ) THEN
	615	IF ( DEBUG ) PRINT *, 'PSFC'
	616	IF (ALLOCATED(psfc)) deallocate(psfc)
	617	allocate (psfc(iweg-1, isng-1, times_in_file ))
	618	IF (ALLOCATED(data1)) deallocate(data1)
	619	allocate (data1(iweg-1, isng-1, 1, times_in_file ))
	620	rcode = nf_inq_varid ( ncid, "PSFC", i )
	621	rcode = nf_get_var_real ( ncid, i, data1 )
	622	psfc(:,:,:) = data1(:,:,1,:)
	623	deallocate (data1)
	624	ENDIF
	625	!
	626	! --> TOPOGRAPHY
	627	!
	628	IF ( LINLOG .ne. 3 ) THEN
	629	IF ( DEBUG ) PRINT *, 'TOPO'
	630	IF (ALLOCATED(ter)) deallocate(ter)
	631	allocate (ter(iweg-1, isng-1))
	632	IF (ALLOCATED(data1)) deallocate(data1)
	633	allocate (data1(iweg-1, isng-1, 1, times_in_file ))
	634	rcode = nf_inq_varid ( ncid, "HGT", i )
	635	rcode = nf_get_var_real ( ncid, i, data1 )
	636	ter(:,:) = data1(:,:,1,1)
	637	IF ( LINLOG .ne. 4 ) deallocate (data1)
	638	ENDIF
	639
	640	! IF (ALLOCATED(qv)) deallocate(qv)
	641	! allocate (qv(iweg-1, isng-1, ibtg-1, times_in_file ))
	642	! rcode = nf_inq_varid ( ncid, "QVAPOR", i )
	643	! rcode = nf_get_var_real ( ncid, i, qv )
	644	!
	645	! NOT SUITABLE for MARS
	646	!
	647	! IF (ALLOCATED(rh)) deallocate(rh)
	648	! allocate (rh(iweg-1, isng-1, ibtg-1, times_in_file ))
	649	! IF (ALLOCATED(data1)) deallocate(data1)
	650	! IF (ALLOCATED(data2)) deallocate(data2)
	651	! allocate (data1(iweg-1, isng-1, ibtg-1, times_in_file ))
	652	! allocate (data2(iweg-1, isng-1, ibtg-1, times_in_file ))
	653	! data1 = 10.0.6112exp(17.67*(tk-273.16)/(TK-29.65))
	654	! data2 = 0.622data1/(0.01 pres_field - (1.-0.622)*data1)
	655	! rh = 100.*AMAX1(AMIN1(qv/data2,1.0),0.0)
	656	! deallocate (data1)
	657	! deallocate (data2)
	658	!
	659	! --> USEFUL for INTERPOLATION
	660	!
	661	!!! at that point pres_field is not used for calculation anymore
	662	!!! it is used to set the initial coordinate for interpolation
	663	IF ( LINLOG == 3) THEN
	664	IF (ALLOCATED(pres_field)) deallocate(pres_field)
	665	allocate ( pres_field(iweg-1, isng-1, ibtg-1, times_in_file) )
	666	pres_field = phb(:,:,1:ibtg-1,:) / grav
	667	ENDIF
	668	IF ( LINLOG == 4) THEN
	669	IF (ALLOCATED(pres_field)) deallocate(pres_field)
	670	allocate ( pres_field(iweg-1, isng-1, ibtg-1, times_in_file))
	671	DO kk = 1, ibtg-1
	672	pres_field(:,:,kk,:) = - data1(:,:,1,:) + phb(:,:,kk,:) / grav
	673	ENDDO
	674	deallocate (data1)
	675	PRINT *, pres_field(10,10,:,1)
	676	ENDIF
	677
	678
	679	IF (ALLOCATED(pres_stagU)) deallocate(pres_stagU)
	680	IF (ALLOCATED(pres_stagV)) deallocate(pres_stagV)
	681	IF ( DEBUG ) PRINT *, 'STAGGERED COORDINATES'
	682	allocate (pres_stagU(iweg, isng-1, ibtg-1, times_in_file ))
	683	allocate (pres_stagV(iweg-1, isng, ibtg-1, times_in_file ))
	684	pres_stagU(1,:,:,:) = pres_field(1,:,:,:)
	685	pres_stagU(iweg,:,:,:) = pres_field(iweg-1,:,:,:)
	686	pres_stagU(2:iweg-1,:,:,:) = (pres_field(1:iweg-2,:,:,:) + pres_field(2:iweg-1,:,:,:))*.5
	687	pres_stagV(:,1,:,:) = pres_field(:,1,:,:)
	688	pres_stagV(:,isng,:,:) = pres_field(:,isng-1,:,:)
	689	pres_stagV(:,2:isng-1,:,:) = (pres_field(:,1:isng-2,:,:) + pres_field(:,2:isng-1,:,:))*.5
	690
	691	!-----------------------------
	692	!-----------------------------
	693	! TRAIN FILE
	694	!-----------------------------
	695	!-----------------------------
	696	IF (debug) THEN
	697	write(6,*)
	698	write(6,*)
	699	write(6,*) "FILE variables:"
	700	ENDIF
	701	jvar = 0
	702	loop_variables : DO ivar = 1, nvars
	703
	704	rcode = nf_inq_var(ncid, ivar, cval, itype, idm, ishape, natt)
	705
	706	!!! Do we want this variable ?
	707	! IF ( trim(cval) == 'P' .OR. trim(cval) == 'PB' ) CYCLE loop_variables
	708	! IF ( trim(cval) == 'PH' .OR. trim(cval) == 'PHB' ) CYCLE loop_variables
	709	!IF ( trim(cval) == 'PTOT' ) CYCLE loop_variables ! PTOT: no
	710	IF ( trim(cval) == 'PHTOT' ) CYCLE loop_variables ! PHTOT: no
	711	IF ( trim(cval) == 'T' ) CYCLE loop_variables ! T: no (tk and tpot calculated above)
	712	IF ( unstagger_grid .AND. (INDEX(cval,'_U') /= 0) ) CYCLE loop_variables !!! no sense in keeping these
	713	IF ( unstagger_grid .AND. (INDEX(cval,'_V') /= 0) ) CYCLE loop_variables !!! no sense in keeping these
	714	IF ( INDEX(process,'all') == 0 ) THEN
	715	!!! Only want some variables - see which
	716	dummy = ","//trim(cval)//","
	717	is_there = INDEX(process_these_fields,trim(dummy))
	718	IF ( is_there == 0 ) THEN
	719	IF ( debug ) print*,"NOTE: ", trim(cval), " - Not requested"
	720	CYCLE loop_variables !!! don't want this one
	721	ENDIF
	722	ENDIF
	723
	724	IF ( idm >= 4 .AND. itype == 4 ) THEN
	725	print*,"NOTE: We cannot deal with 3D integers - maybe later"
	726	CYCLE loop_variables
	727	ENDIF
	728	IF ( itype == 6 ) THEN
	729	print*,"NOTE: We cannot deal with double precision data - maybe later"
	730	CYCLE loop_variables
	731	ENDIF
	732	IF ( itype == 2 .OR. itype == 4 .OR. itype == 5 ) THEN
	733	!!! OK I know what to do this this
	734	ELSE
	735	print*,"NOTE: Do not understand this data type ", itype, " skip field."
	736	CYCLE loop_variables
	737	ENDIF
	738
	739	!IF ( trim(cval) == 'U' ) cval = 'UU'
	740	!IF ( trim(cval) == 'V' ) cval = 'VV'
	741	!IF ( trim(cval) == 'T' ) cval = 'TT'
	742
	743	!!! OK - we want this - lets continue
	744	jvar = jvar + 1
	745	jshape = 0
	746	interpolate = .FALSE.
	747	fix_meta_stag = .FALSE.
	748	rcode = nf_redef(mcid)
	749	DO ii = 1, idm
	750	test_dim_name = dnamei(ishape(ii))
	751	IF ( test_dim_name == 'bottom_top' .OR. test_dim_name == 'bottom_top_stag' ) THEN
	752	IF ( test_dim_name == 'bottom_top_stag' ) fix_meta_stag = .TRUE.
	753	IF (LINLOG .le. 2) test_dim_name = 'pressure'
	754	IF (LINLOG .eq. 3) test_dim_name = 'altitude'
	755	IF (LINLOG .eq. 4) test_dim_name = 'altitude_abg'
	756	interpolate = .TRUE.
	757	ENDIF
	758	IF ( unstagger_grid .AND. test_dim_name == 'west_east_stag' ) THEN
	759	test_dim_name = 'west_east'
	760	fix_meta_stag = .TRUE.
	761	ENDIF
	762	IF ( unstagger_grid .AND. test_dim_name == 'south_north_stag' ) THEN
	763	test_dim_name = 'south_north'
	764	fix_meta_stag = .TRUE.
	765	ENDIF
	766	DO jj = 1,j
	767	IF ( test_dim_name == dnamej(jj) ) THEN
	768	jshape(ii) = jj
	769	ENDIF
	770	ENDDO
	771
	772	IF ( jshape(ii) == 0 ) THEN
	773	j = j + 1
	774	jshape(ii) = j
	775	dnamej(j) = dnamei(ishape(ii))
	776	dvalj(j) = dvali(ishape(ii))
	777	rcode = nf_def_dim(mcid, dnamej(j), dvalj(j), j)
	778	ENDIF
	779	ENDDO
	780	rcode = nf_def_var(mcid, cval, itype, idm, jshape, jvar)
	781	!print *, 'mcid', mcid
	782	!print *, 'cval', cval
	783	!print *, 'itype', itype
	784	!print *, 'idm', idm
	785	!print *, 'jshape', jshape
	786	!print *, 'jvar', jvar
	787
	788	DO na = 1, natt
	789	rcode = nf_inq_attname(ncid, ivar, na, cname)
	790	IF ( fix_meta_stag .AND. trim(cname) == 'stagger' ) THEN
	791	att_text = "-"
	792	ilen = len_trim(att_text)
	793	rcode = nf_put_att_text(mcid, jvar, cname, ilen, att_text(1:ilen) )
	794	ELSEIF ( fix_meta_stag .AND. trim(cname) == 'coordinates' ) THEN
	795	att_text = "XLONG XLAT"
	796	ilen = len_trim(att_text)
	797	rcode = nf_put_att_text(mcid, jvar, cname, ilen, att_text(1:ilen) )
	798	ELSE
	799	rcode = nf_copy_att(ncid, ivar, cname, mcid, jvar)
	800	ENDIF
	801	ENDDO
	802
	803	IF ( extrapolate == 0 ) THEN
	804	rcode = nf_put_att_real(mcid, jvar, 'missing_value', NF_FLOAT, 1, MISSING )
	805	ENDIF
	806
	807	rcode = nf_enddef(mcid)
	808
	809	!
	810	! GET THE DIMS FOR INPUT AND OUTPUT FROM THE SHAPE
	811	!
	812	dims_in = 1
	813	dims_out = 1
	814	DO ii = 1,idm
	815	dims_in(ii) = dvali(ishape(ii))
	816	dims_out(ii) = dvalj(jshape(ii))
	817	ENDDO
	818	IF (debug) write(6,*) 'VAR: ',trim(cval)
	819	IF (debug) THEN
	820	write(6,*) ' DIMS IN: ',dims_in
	821	write(6,*) ' DIMS OUT: ',dims_out
	822	ENDIF
	823
	824	!
	825	! ALLOCATE THE INPUT AND OUTPUT ARRAYS
	826	! READ THE DATA FROM INPUT FILE
	827	!
	828	IF (itype == 2) THEN ! character
	829	allocate (text(dims_in(1), dims_in(2), dims_in(3), dims_in(4)))
	830	rcode = nf_get_var_text(ncid, ivar, text)
	831	rcode = nf_put_vara_text (mcid, jvar, start_dims, dims_in, text)
	832	IF (debug) write(6,*) ' SAMPLE VALUE = ',text(:,1,1,1)
	833	deallocate (text)
	834
	835	ELSEIF (itype == 4) THEN ! integer
	836	allocate (idata1(dims_in(1), dims_in(2), dims_in(3), dims_in(4)))
	837	rcode = nf_get_var_int(ncid, ivar, idata1)
	838	rcode = nf_put_vara_int (mcid, jvar, start_dims, dims_in, idata1)
	839	IF (debug) write(6,*) ' SAMPLE VALUE = ',idata1(dims_in(1)/2,dims_in(2)/2,1,1)
	840	deallocate (idata1)
	841
	842	ELSEIF (itype == 5) THEN ! real
	843	allocate (data1(dims_in(1), dims_in(2), dims_in(3), dims_in(4)))
	844	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
	845	rcode = nf_get_var_real(ncid, ivar, data1)
	846
	847	IF (idm >= 4 .AND. interpolate) THEN
	848	IF (debug) write(6,*) ' THIS IS A FIELD WE NEED TO INTERPOLATE'
	849
	850	IF ( dims_in(1) == iweg .AND. .not. unstagger_grid ) THEN
	851	CALL interp (data2, data1, pres_stagU, interp_levels, psfc, ter, tk, qv, &
	852	iweg, isng-1, ibtg-1, dims_in(4), &
	853	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	854	ELSEIF ( dims_in(2) == isng .AND. .not. unstagger_grid ) THEN
	855	CALL interp (data2, data1, pres_stagV, interp_levels, psfc, ter, tk, qv, &
	856	iweg-1, isng, ibtg-1, dims_in(4), &
	857	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	858	ELSEIF ( dims_in(1) == iweg .AND. unstagger_grid ) THEN
	859	allocate (data3(iweg-1, isng-1, ibtg-1, dims_in(4)))
	860	data3(1:iweg-1,:,:,:) = (data1(1:iweg-1,:,:,:) + data1(2:iweg,:,:,:)) * .5
	861	CALL interp (data2, data3, pres_field, interp_levels, psfc, ter, tk, qv, &
	862	iweg-1, isng-1, ibtg-1, dims_in(4), &
	863	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	864	deallocate(data3)
	865	ELSEIF ( dims_in(2) == isng .AND. unstagger_grid ) THEN
	866	allocate (data3(iweg-1, isng-1, ibtg-1, dims_in(4)))
	867	data3(:,1:isng-1,:,:) = (data1(:,1:isng-1,:,:) + data1(:,2:isng,:,:)) * .5
	868	CALL interp (data2, data3, pres_field, interp_levels, psfc, ter, tk, qv, &
	869	iweg-1, isng-1, ibtg-1, dims_in(4), &
	870	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	871	deallocate(data3)
	872	ELSEIF ( dims_in(3) == ibtg ) THEN
	873	allocate (data3(iweg-1, isng-1, ibtg-1, dims_in(4)))
	874	IF (debug) PRINT *, 'VERTICAL UNSTAGGERING'
	875	data3(:,:,1:ibtg-1,:) = (data1(:,:,1:ibtg-1,:) + data1(:,:,2:ibtg,:)) * .5
	876	CALL interp (data2, data3, pres_field, interp_levels, psfc, ter, tk, qv, &
	877	dims_in(1), dims_in(2), ibtg-1, dims_in(4), &
	878	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	879	deallocate(data3)
	880	ELSE
	881	CALL interp (data2, data1, pres_field, interp_levels, psfc, ter, tk, qv, &
	882	dims_in(1), dims_in(2), dims_in(3), dims_in(4), &
	883	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	884	END IF
	885
	886	IF (debug) write(6,*) ' SAMPLE VALUE IN = ',data1(dims_in(1)/2,dims_in(2)/2,1,1)
	887	IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',data2(dims_out(1)/2,dims_out(2)/2,1,1)
	888
	889	ELSEIF (idm == 3 .AND. unstagger_grid ) THEN
	890	IF ( dims_in(1) == iweg ) THEN
	891	data2(1:iweg-1,:,:,:) = (data1(1:iweg-1,:,:,:) + data1(2:iweg,:,:,:)) * .5
	892	ELSEIF ( dims_in(2) == isng ) THEN
	893	data2(:,1:isng-1,:,:) = (data1(:,1:isng-1,:,:) + data1(:,2:isng,:,:)) * .5
	894	ELSE
	895	data2 = data1
	896	ENDIF
	897	IF (debug) write(6,*) ' SAMPLE VALUE = ',data1(dims_in(1)/2,dims_in(2)/2,1,1)
	898
	899	ELSE
	900	data2 = data1
	901	IF (debug) write(6,*) ' SAMPLE VALUE = ',data1(dims_in(1)/2,dims_in(2)/2,1,1)
	902
	903	ENDIF
	904
	905	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
	906	!print *, 'mcid ', mcid
	907	!print *, 'jvar ', jvar
	908	!print *, 'start_dims ', start_dims
	909	!print *, 'dims_out ', dims_out
	910
	911	deallocate (data1)
	912	deallocate (data2)
	913
	914	ENDIF
	915
	916	ENDDO loop_variables
	917
	918	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	919	!!!!!! ADDITIONAL VARIABLES !!!!!!
	920	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	921	IF ( debug ) print*," "
	922	IF ( debug ) print*,"Calculating some diagnostics"
	923	interpolate = .TRUE.
	924
	925	!!! NB: what follows is not useful because we'd like diagnostics for each history timestep
	926	jshape = 0
	927	DO ii = 1, 4
	928	IF ( ii == 1 ) test_dim_name = 'west_east'
	929	IF ( ii == 2 ) test_dim_name = 'south_north'
	930	IF (( ii == 3 ) .and. (LINLOG .le. 2)) test_dim_name = 'pressure'
	931	IF (( ii == 3 ) .and. (LINLOG .eq. 3)) test_dim_name = 'altitude'
	932	IF (( ii == 3 ) .and. (LINLOG .eq. 4)) test_dim_name = 'altitude_abg'
	933	IF ( ii == 4 ) test_dim_name = 'Time'
	934	DO jj = 1,j
	935	IF ( test_dim_name == dnamej(jj) ) THEN
	936	jshape(ii) = jj
	937	ENDIF
	938	ENDDO
	939	IF ( jshape(ii) == 0 ) THEN
	940	j = j + 1
	941	jshape(ii) = j
	942	dnamej(j) = dnamei(ishape(ii))
	943	dvalj(j) = dvali(ishape(ii))
	944	rcode = nf_def_dim(mcid, dnamej(j), dvalj(j), j)
	945	ENDIF
	946	ENDDO
	947	dims_in = 1
	948	dims_out = 1
	949	DO ii = 1,4
	950	dims_out(ii) = dvalj(jshape(ii))
	951	print *, dims_out(ii)
	952	ENDDO
	953	!!! NB: what follows is useful because we'd like diagnostics for each history timestep
	954	dims_in = dims_out
	955
	956	!IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'PRES') /= 0 ) THEN
	957	! !!! PRES
	958	! jvar = jvar + 1
	959	! CALL def_var (mcid, jvar, "PRES", 5, 4, jshape, "XZY", "Pressure ", "Pa", "-", "XLONG XLAT")
	960	! IF (debug) THEN
	961	! write(6,*) 'VAR: PRES'
	962	! write(6,*) ' DIMS OUT: ',dims_out
	963	! ENDIF
	964	! rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, pres_out)
	965	! IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',pres_out(dims_out(1)/2,dims_out(2)/2,1,1)
	966	!ENDIF
	967
	968	!
	969	! EN FAIT IL FAUDRAIT QUE LE jshape SOIT BIEN DEFINI POUR CE QUI SUIT
	970	!
	971
	972	!
	973	! OUTPUT REGULAR TEMPERATURE
	974	!
	975	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'tk') /= 0 ) THEN
	976	jvar = jvar + 1
	977	CALL def_var (mcid, jvar, "tk ", 5, 4, jshape, "XZY", "Temperature ", "K ", "-", "XLONG XLAT", MISSING)
	978	IF (debug) THEN
	979	write(6,*) 'VAR: tk'
	980	write(6,*) ' DIMS OUT: ',dims_out
	981	ENDIF
	982	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
	983	CALL interp (data2, tk, pres_field, interp_levels, psfc, ter, tk, qv, &
	984	iweg-1, isng-1, ibtg-1, dims_in(4), &
	985	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	986	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
	987	IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',data2(dims_out(1)/2,dims_out(2)/2,1,1)
	988	deallocate(data2)
	989	ENDIF
	990
	991	!
	992	! OUTPUT ROTATED WINDS
	993	!
	994	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'uvmet') /= 0 ) THEN
	995	jvar = jvar + 1
	996	CALL def_var (mcid, jvar, "Um ", 5, 4, jshape, "XZY", "U rotated wind ", "K ", "-", "XLONG XLAT", MISSING)
	997	IF (debug) THEN
	998	write(6,*) 'VAR: u (rotated)'
	999	write(6,*) ' DIMS OUT: ',dims_out
	1000	ENDIF
	1001	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
	1002	CALL interp (data2, umet, pres_field, interp_levels, psfc, ter, tk, qv, &
	1003	iweg-1, isng-1, ibtg-1, dims_in(4), &
	1004	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	1005	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
	1006	IF (debug) write(6,*) ' SAMPLE VALUE OUT=',data2(dims_out(1)/2,dims_out(2)/2,1,1)
	1007	deallocate(data2)
	1008
	1009	jvar = jvar + 1
	1010	CALL def_var (mcid, jvar, "Vm ", 5, 4, jshape, "XZY", "V rotated wind ", "K ", "-", "XLONG XLAT", MISSING)
	1011	IF (debug) THEN
	1012	write(6,*) 'VAR: v (rotated)'
	1013	write(6,*) ' DIMS OUT: ',dims_out
	1014	ENDIF
	1015	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
	1016	CALL interp (data2, vmet, pres_field, interp_levels, psfc, ter, tk, qv, &
	1017	iweg-1, isng-1, ibtg-1, dims_in(4), &
	1018	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	1019	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
	1020	IF (debug) write(6,*) ' SAMPLE VALUE OUT=',data2(dims_out(1)/2,dims_out(2)/2,1,1)
	1021	deallocate(data2)
	1022	ENDIF
	1023
	1024	!
	1025	! OUTPUT POTENTIAL TEMPERATURE
	1026	!
	1027	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'tpot') /= 0 ) THEN
	1028	jvar = jvar + 1
	1029	CALL def_var (mcid, jvar, "tpot", 5, 4, jshape, "XZY", "Potential Temperature ", "K ", "-", "XLONG XLAT", MISSING)
	1030	IF (debug) THEN
	1031	write(6,*) 'VAR: tpot'
	1032	write(6,*) ' DIMS OUT: ',dims_out
	1033	ENDIF
	1034	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
	1035	CALL interp (data2, tpot, pres_field, interp_levels, psfc, ter, tk, qv, &
	1036	iweg-1, isng-1, ibtg-1, dims_in(4), &
	1037	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	1038	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
	1039	IF (debug) write(6,*) ' SAMPLE VALUE OUT =',data2(dims_out(1)/2,dims_out(2)/2,1,1)
	1040	deallocate(data2)
	1041	ENDIF
	1042
	1043	!
	1044	! OUTPUT GEOPOTENTIAL HEIGHT
	1045	!
	1046	IF (LINLOG .le. 2) THEN
	1047	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'GHT') /= 0 ) THEN
	1048	jvar = jvar + 1
	1049	CALL def_var (mcid, jvar, "GHT ", 5, 4, jshape, "XZY", "Geopotential Height", "m ", "-", "XLONG XLAT", MISSING)
	1050	IF (debug) THEN
	1051	write(6,*) 'VAR: GHT'
	1052	write(6,*) ' DIMS OUT: ',dims_out
	1053	ENDIF
	1054	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
	1055	CALL interp (data2, phb(:,:,1:ibtg-1,:), pres_field, interp_levels, psfc, ter, tk, qv, &
	1056	iweg-1, isng-1, ibtg-1, dims_in(4), &
	1057	num_metgrid_levels, LINLOG, extrapolate, .TRUE., MISSING)
	1058	data2 = data2/grav
	1059	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
	1060	IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',data2(dims_out(1)/2,dims_out(2)/2,1,1)
	1061	deallocate(data2)
	1062	ENDIF
	1063	ELSE
	1064	PRINT *, 'Geopotential height is actually vertical coordinate'
	1065	ENDIF
	1066
	1067	! IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'RH') /= 0 ) THEN
	1068	! !!! RH
	1069	! jvar = jvar + 1
	1070	! CALL def_var (mcid, jvar, "RH ", 5, 4, jshape, "XZY", "Relative Humidity ", "% ", "-", "XLONG XLAT")
	1071	! IF (debug) THEN
	1072	! write(6,*) 'VAR: RH'
	1073	! write(6,*) ' DIMS OUT: ',dims_out
	1074	! ENDIF
	1075	! allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
	1076	! CALL interp (data2, rh, pres_field, interp_levels, psfc, ter, tk, qv, &
	1077	! iweg-1, isng-1, ibtg-1, dims_in(4), &
	1078	! num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
	1079	! WHERE ( rh < 0.0 )
	1080	! rh = 0.0
	1081	! ENDWHERE
	1082	! WHERE ( rh > 100.0 )
	1083	! rh = 100.0
	1084	! ENDWHERE
	1085	! rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
	1086	! IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',data2(dims_out(1)/2,dims_out(2)/2,1,1)
	1087	! deallocate(data2)
	1088	! ENDIF
	1089
	1090
	1091	!=====================================================================================
	1092	!=====================================================================================
	1093	!
	1094	! VERTICAL COORDINATES
	1095	!
	1096	jvar = jvar + 1
	1097	jshape(:)=0
	1098	jshape(1)=2
	1099	CALL def_var (mcid, jvar, "vert", 5, 1, jshape, " Z ", "Vert. coord. ","m ", "-", "XLONG XLAT", MISSING)
	1100	start_dims(1)=1
	1101	start_dims(2)=1
	1102	start_dims(3)=1
	1103	start_dims(4)=1
	1104	dims_out(1)=dims_out(3)
	1105	dims_out(2)=1
	1106	dims_out(3)=1
	1107	dims_out(4)=1
	1108	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
	1109	do kk=1,dims_out(1)
	1110	data2(kk,1,1,1) = interp_levels(kk)
	1111	enddo
	1112	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
	1113	deallocate(data2)
	1114	!=====================================================================================
	1115	!=====================================================================================
	1116
	1117
	1118
	1119	!
	1120	! CLOSE FILES ON EXIT
	1121	!
	1122	rcode = nf_close(ncid)
	1123	rcode = nf_close(mcid)
	1124	write(6,*)
	1125	ENDDO
	1126
	1127	!
	1128	! WELL...
	1129	!
	1130	write(6,*)
	1131	write(6,*) "##########################################"
	1132	write(6,*) " END of API PROGRAM - SUCCESS so far"
	1133	write(6,*) "##########################################"
	1134
	1135	END PROGRAM api
	1136	!---------------------------------------------------------------------
	1137	!---------------------------------------------------------------------
	1138	!---------------------------------------------------------------------
	1139	SUBROUTINE handle_err(rcode)
	1140	INTEGER rcode
	1141	write(6,*) 'Error number ',rcode
	1142	stop
	1143	END SUBROUTINE
	1144	!---------------------------------------------------------------------
	1145	SUBROUTINE interp (data_out, data_in, pres_field, interp_levels, psfc, ter, tk, qv, ix, iy, iz, it, &
	1146	num_metgrid_levels, LINLOG, extrapolate, GEOPT, MISSING)
	1147
	1148	INTEGER :: ix, iy, iz, it
	1149	INTEGER :: num_metgrid_levels, LINLOG
	1150	REAL, DIMENSION(ix, iy, num_metgrid_levels, it) :: data_out
	1151	REAL, DIMENSION(ix, iy, iz, it) :: data_in, pres_field, tk, qv
	1152	REAL, DIMENSION(ix, iy, it) :: psfc
	1153	REAL, DIMENSION(ix, iy) :: ter
	1154	REAL, DIMENSION(num_metgrid_levels) :: interp_levels
	1155
	1156	INTEGER :: i, j, itt, k, kk, kin
	1157	REAL, DIMENSION(num_metgrid_levels) :: data_out1D
	1158	REAL, DIMENSION(iz) :: data_in1D, pres_field1D
	1159	INTEGER :: extrapolate
	1160	REAL :: MISSING
	1161	REAL, DIMENSION(ix, iy, num_metgrid_levels, it) :: N
	1162	REAL :: sumA, sumN, AVE_geopt
	1163	LOGICAL :: GEOPT
	1164
	1165	N = 1.0
	1166
	1167	do itt = 1, it
	1168	!PRINT *, 'TIME... ', itt, ' in ', iz, ' out ', num_metgrid_levels
	1169	do j = 1, iy
	1170	do i = 1, ix
	1171	data_in1D(:) = data_in(i,j,:,itt)
	1172	pres_field1D(:) = pres_field(i,j,:,itt)
	1173	IF ( LINLOG .le. 2 ) THEN
	1174	CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING)
	1175	ELSE
	1176	CALL interp_1d (data_in1D, pres_field1D, iz, data_out1D, interp_levels, num_metgrid_levels, 'z', MISSING)
	1177	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1178	!interp_1d (data_in1D, pres_field1D, num_metgrid_levels, data_out1D, interp_levels, iz, 'z')
	1179	!CALL interp_1d( data_in_1d, z_data_1d, bottom_top_dim, data_out_1d, z_levs, number_of_zlevs, vertical_type)
	1180	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1181	ENDIF
	1182	data_out(i,j,:,itt) = data_out1D(:)
	1183	end do
	1184	end do
	1185	end do
	1186	!PRINT *, 'ok'
	1187
	1188	! Fill in missing values
	1189	IF ( extrapolate == 0 ) RETURN !! no extrapolation - we are out of here
	1190
	1191
	1192
	1193	IF ( LINLOG .ge. 0 ) RETURN !! TEMPORARY: no extrapolation
	1194	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1195	!! CAUTION BELOW: METHOD NOT ADAPTED TO MARS
	1196	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1197
	1198	!!! MARS MARS MARS
	1199	!expon=287.04*.0065/9.81
	1200	expon=192.*.0045/3.72
	1201
	1202	! First find where about 400 hPa is located
	1203	kk = 0
	1204	find_kk : do k = 1, num_metgrid_levels
	1205	kk = k
	1206	if ( interp_levels(k) <= 40000. ) exit find_kk
	1207	end do find_kk
	1208
	1209
	1210	IF ( GEOPT ) THEN !! geopt is treated different below ground
	1211
	1212	do itt = 1, it
	1213	do k = 1, kk
	1214	do j = 1, iy
	1215	do i = 1, ix
	1216	IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN
	1217
	1218	! We are below the first model level, but above the ground
	1219
	1220	!!! MARS MARS
	1221	data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))ter(i,j)3.72 + &
	1222	(psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) / &
	1223	(psfc(i,j,itt) - pres_field(i,j,1,itt))
	1224
	1225	ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN
	1226
	1227	! We are below both the ground and the lowest data level.
	1228
	1229	! First, find the model level that is closest to a "target" pressure
	1230	! level, where the "target" pressure is delta-p less that the local
	1231	! value of a horizontally smoothed surface pressure field. We use
	1232	! delta-p = 150 hPa here. A standard lapse rate temperature profile
	1233	! passing through the temperature at this model level will be used
	1234	! to define the temperature profile below ground. This is similar
	1235	! to the Benjamin and Miller (1990) method, except that for
	1236	! simplicity, they used 700 hPa everywhere for the "target" pressure.
	1237	! Code similar to what is implemented in RIP4
	1238
	1239
	1240	!!! MARS MARS MARS
	1241	ptarget = (psfc(i,j,itt)*.01) - 1.50 !150.
	1242	dpmin=1.e2 !1.e4
	1243	kupper = 0
	1244	loop_kIN : do kin=iz,1,-1
	1245	kupper = kin
	1246	dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget )
	1247	if (dp.gt.dpmin) exit loop_kIN
	1248	dpmin=min(dpmin,dp)
	1249	enddo loop_kIN
	1250
	1251	pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt))
	1252	zbot=min(data_in(i,j,1,itt)/3.72,ter(i,j)) !!MARS MARS
	1253
	1254	tbotextrap=tk(i,j,kupper,itt)(pbot/pres_field(i,j,kupper,itt))*expon
	1255	tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
	1256
	1257	data_out(i,j,k,itt) = (zbot+tvbotextrap/.0045(1.-(interp_levels(k)/pbot)expon))3.72 !!MARS MARS
	1258
	1259	ENDIF
	1260	enddo
	1261	enddo
	1262	enddo
	1263	enddo
	1264
	1265
	1266	!!! Code for filling missing data with an average - we don't want to do this
	1267	!!do itt = 1, it
	1268	!!loop_levels : do k = 1, num_metgrid_levels
	1269	!!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
	1270	!!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
	1271	!!IF ( sumN == 0. ) CYCLE loop_levels
	1272	!!AVE_geopt = sumA/sumN
	1273	!!WHERE ( data_out(:,:,k,itt) == MISSING )
	1274	!!data_out(:,:,k,itt) = AVE_geopt
	1275	!!END WHERE
	1276	!!end do loop_levels
	1277	!!end do
	1278
	1279	END IF
	1280
	1281	!!! All other fields and geopt at higher levels come here
	1282	do itt = 1, it
	1283	do j = 1, iy
	1284	do i = 1, ix
	1285	do k = 1, kk
	1286	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt)
	1287	end do
	1288	do k = kk+1, num_metgrid_levels
	1289	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt)
	1290	end do
	1291	end do
	1292	end do
	1293	end do
	1294
	1295
	1296	END SUBROUTINE interp
	1297	!------------------------------------------------------------------------------
	1298	!--------------------------------------------------------
	1299	SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING)
	1300
	1301	! Modified from int2p - NCL code
	1302	! routine to interpolate from one set of pressure levels
	1303	! . to another set using linear or ln(p) interpolation
	1304	!
	1305	! NCL: xout = int2p (pin,xin,pout,linlog)
	1306	! This code was originally written for a specific purpose.
	1307	! . Several features were added for incorporation into NCL's
	1308	! . function suite including linear extrapolation.
	1309	!
	1310	! nomenclature:
	1311	!
	1312	! . ppin - input pressure levels. The pin can be
	1313	! . be in ascending or descending order
	1314	! . xxin - data at corresponding input pressure levels
	1315	! . npin - number of input pressure levels >= 2
	1316	! . ppout - output pressure levels (input by user)
	1317	! . same (ascending or descending) order as pin
	1318	! . xxout - data at corresponding output pressure levels
	1319	! . npout - number of output pressure levels
	1320	! . linlog - if abs(linlog)=1 use linear interp in pressure
	1321	! . if abs(linlog)=2 linear interp in ln(pressure)
	1322	! . missing- missing data code.
	1323
	1324	! ! input types
	1325	INTEGER :: npin,npout,linlog,ier
	1326	real :: ppin(npin),xxin(npin),ppout(npout)
	1327	real :: MISSING
	1328	logical :: AVERAGE
	1329	! ! output
	1330	real :: xxout(npout)
	1331	INTEGER :: j1,np,nl,nin,nlmax,nplvl
	1332	INTEGER :: nlsave,np1,no1,n1,n2,nlstrt
	1333	real :: slope,pa,pb,pc
	1334
	1335	! automatic arrays
	1336	real :: pin(npin),xin(npin),p(npin),x(npin)
	1337	real :: pout(npout),xout(npout)
	1338
	1339
	1340	xxout = MISSING
	1341	pout = ppout
	1342	p = ppin
	1343	x = xxin
	1344	nlmax = npin
	1345
	1346	! exact p-level matches
	1347	nlstrt = 1
	1348	nlsave = 1
	1349	do np = 1,npout
	1350	xout(np) = MISSING
	1351	do nl = nlstrt,nlmax
	1352	if (pout(np).eq.p(nl)) then
	1353	xout(np) = x(nl)
	1354	nlsave = nl + 1
	1355	go to 10
	1356	end if
	1357	end do
	1358	10 nlstrt = nlsave
	1359	end do
	1360
	1361	if (LINLOG.eq.1) then
	1362	do np = 1,npout
	1363	do nl = 1,nlmax - 1
	1364	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
	1365	slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
	1366	xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
	1367	end if
	1368	end do
	1369	end do
	1370	elseif (LINLOG.eq.2) then
	1371	do np = 1,npout
	1372	do nl = 1,nlmax - 1
	1373	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
	1374	pa = log(p(nl))
	1375	pb = log(pout(np))
	1376	! special case: in case someone inadvertently enter p=0.
	1377	if (p(nl+1).gt.0.d0) then
	1378	pc = log(p(nl+1))
	1379	else
	1380	pc = log(1.d-4)
	1381	end if
	1382
	1383	slope = (x(nl)-x(nl+1))/ (pa-pc)
	1384	xout(np) = x(nl+1) + slope* (pb-pc)
	1385	end if
	1386	end do
	1387	end do
	1388	end if
	1389
	1390
	1391	! place results in the return array;
	1392	xxout = xout
	1393
	1394	END SUBROUTINE int1D
	1395
	1396	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1397	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1398	SUBROUTINE interp_1d( a, xa, na, b, xb, nb, vertical_type, MISSING)
	1399
	1400	implicit none
	1401
	1402	! Arguments
	1403	integer, intent(in) :: na, nb
	1404	real, intent(in), dimension(na) :: a, xa
	1405	real, intent(in), dimension(nb) :: xb
	1406	real, intent(out), dimension(nb) :: b
	1407	character (len=1) :: vertical_type
	1408
	1409	!! Local variables
	1410	real :: MISSING
	1411	integer :: n_in, n_out
	1412	real :: w1, w2
	1413	logical :: interp
	1414
	1415	IF ( vertical_type == 'p' ) THEN
	1416
	1417	DO n_out = 1, nb
	1418
	1419	b(n_out) = MISSING
	1420	interp = .false.
	1421	n_in = 1
	1422
	1423	DO WHILE ( (.not.interp) .and. (n_in < na) )
	1424	IF( (xa(n_in) >= xb(n_out)) .and. &
	1425	(xa(n_in+1) <= xb(n_out)) ) THEN
	1426	interp = .true.
	1427	w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
	1428	w2 = 1. - w1
	1429	b(n_out) = w1a(n_in) + w2a(n_in+1)
	1430	END IF
	1431	n_in = n_in +1
	1432	ENDDO
	1433
	1434	ENDDO
	1435
	1436	ELSE
	1437
	1438	DO n_out = 1, nb
	1439
	1440	b(n_out) = MISSING
	1441	interp = .false.
	1442	n_in = 1
	1443
	1444	DO WHILE ( (.not.interp) .and. (n_in < na) )
	1445	IF( (xa(n_in) <= xb(n_out)) .and. &
	1446	(xa(n_in+1) >= xb(n_out)) ) THEN
	1447	interp = .true.
	1448	w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
	1449	w2 = 1. - w1
	1450	b(n_out) = w1a(n_in) + w2a(n_in+1)
	1451	END IF
	1452	n_in = n_in +1
	1453	ENDDO
	1454
	1455	ENDDO
	1456
	1457	END IF
	1458
	1459	END SUBROUTINE interp_1d
	1460	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1461	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1462
	1463	!------------------------------------------------------------------------------
	1464	FUNCTION virtual (tmp,rmix)
	1465	! This function returns virtual temperature in K, given temperature
	1466	! in K and mixing ratio in kg/kg.
	1467
	1468	real :: tmp, rmix, virtual
	1469
	1470	virtual=tmp(0.622+rmix)/(0.622(1.+rmix))
	1471
	1472	END FUNCTION virtual
	1473
	1474	!------------------------------------------------------------------------------
	1475	SUBROUTINE all_spaces ( command , length_of_char )
	1476
	1477	IMPLICIT NONE
	1478
	1479	INTEGER :: length_of_char
	1480	CHARACTER (LEN=length_of_char) :: command
	1481	INTEGER :: loop
	1482
	1483	DO loop = 1 , length_of_char
	1484	command(loop:loop) = ' '
	1485	END DO
	1486
	1487	END SUBROUTINE all_spaces
	1488
	1489	!------------------------------------------------------------------------------
	1490	SUBROUTINE def_var (mcid, jvar, cval, itype, idm, jshape, order, desc, units, stag, coord, missing )
	1491
	1492	IMPLICIT NONE
	1493
	1494	INCLUDE 'netcdf.inc'
	1495
	1496	INTEGER :: mcid, jvar
	1497	CHARACTER (LEN = 4) :: cval
	1498	INTEGER :: itype, idm
	1499	REAL, DIMENSION(6) :: jshape
	1500	CHARACTER (LEN = 3) :: order
	1501	CHARACTER (LEN = 19) :: desc
	1502	CHARACTER (LEN = 2) :: units
	1503	CHARACTER (LEN = 1) :: stag
	1504	CHARACTER (LEN = 10) :: coord
	1505
	1506	INTEGER :: rcode, ilen
	1507	CHARACTER (LEN=30) :: att_text
	1508
	1509	REAL :: missing
	1510
	1511
	1512	IF ( itype == 5 ) THEN
	1513	rcode = nf_redef(mcid)
	1514	rcode = nf_def_var(mcid, trim(cval), NF_REAL, idm, jshape, jvar)
	1515	rcode = nf_put_att_int(mcid, jvar, "FieldType", NF_INT, 1, 104)
	1516	ENDIF
	1517
	1518	att_text = order
	1519	ilen = len_trim(att_text)
	1520	rcode = nf_put_att_text(mcid, jvar, "MemoryOrder", ilen, att_text(1:ilen) )
	1521
	1522	att_text = desc
	1523	ilen = len_trim(att_text)
	1524	rcode = nf_put_att_text(mcid, jvar, "description", ilen, att_text(1:ilen) )
	1525
	1526	att_text = units
	1527	ilen = len_trim(att_text)
	1528	rcode = nf_put_att_text(mcid, jvar, "units", ilen, att_text(1:ilen) )
	1529
	1530	att_text = stag
	1531	ilen = len_trim(att_text)
	1532	rcode = nf_put_att_text(mcid, jvar, "stagger", ilen, att_text(1:ilen) )
	1533
	1534	att_text = coord
	1535	ilen = len_trim(att_text)
	1536	rcode = nf_put_att_text(mcid, jvar, "coordinates", ilen, att_text(1:ilen) )
	1537
	1538	rcode = nf_put_att_real(mcid, jvar, "missing_value", NF_FLOAT, 1, MISSING )
	1539
	1540	rcode = nf_enddef(mcid)
	1541
	1542	END SUBROUTINE def_var
	1543	!------------------------------------------------------------------------------
	1544	!------------------------------------------------------------------------------
	1545	!! Diagnostics: U & V on earth coordinates
	1546	SUBROUTINE calc_uvmet( UUU, VVV, &
	1547	UUUmet, VVVmet, &
	1548	truelat1, truelat2, &
	1549	stand_lon, map_proj, &
	1550	longi, lati, &
	1551	west_east_dim, south_north_dim, bottom_top_dim)
	1552
	1553	IMPLICIT NONE
	1554
	1555	!Arguments
	1556	integer :: west_east_dim, south_north_dim, bottom_top_dim
	1557	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: UUU
	1558	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: VVV
	1559	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: UUUmet
	1560	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: VVVmet
	1561	real, dimension(west_east_dim,south_north_dim) :: longi, lati
	1562	real :: truelat1, truelat2, stand_lon
	1563	integer :: map_proj
	1564
	1565	!Local
	1566	integer :: i, j, k
	1567	real :: cone
	1568	real, dimension(west_east_dim,south_north_dim) :: diff, alpha
	1569
	1570	real, parameter :: PI = 3.141592653589793
	1571	real, parameter :: DEG_PER_RAD = 180./PI
	1572	real, parameter :: RAD_PER_DEG = PI/180.
	1573
	1574	! print *, 'map ', map_proj
	1575	IF ( map_proj .ge. 3 ) THEN ! No need to rotate
	1576	!PRINT *, 'NO NEED TO ROTATE !!!! equivalent to output U,V with unstagger_grid'
	1577	UUUmet(:,:,:) = UUU
	1578	VVVmet(:,:,:) = VVV
	1579	ELSE
	1580	!END IF
	1581
	1582	cone = 1. ! PS
	1583	IF ( map_proj .eq. 1) THEN ! Lambert Conformal mapping
	1584	IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
	1585	cone=(ALOG(COS(truelat1*RAD_PER_DEG))- &
	1586	ALOG(COS(truelat2*RAD_PER_DEG))) / &
	1587	(ALOG(TAN((90.-ABS(truelat1))RAD_PER_DEG0.5 ))- &
	1588	ALOG(TAN((90.-ABS(truelat2))RAD_PER_DEG0.5 )) )
	1589	ELSE
	1590	cone = SIN(ABS(truelat1)*RAD_PER_DEG )
	1591	ENDIF
	1592	END IF
	1593
	1594
	1595	diff = longi - stand_lon
	1596	DO i = 1, west_east_dim
	1597	DO j = 1, south_north_dim
	1598	IF ( diff(i,j) .gt. 180. ) THEN
	1599	diff(i,j) = diff(i,j) - 360.
	1600	END IF
	1601	IF ( diff(i,j) .lt. -180. ) THEN
	1602	diff(i,j) = diff(i,j) + 360.
	1603	END IF
	1604	END DO
	1605	END DO
	1606
	1607	!print *, longi(10,10)
	1608	!print *, lati(10,10)
	1609
	1610
	1611	DO i = 1, west_east_dim
	1612	DO j = 1, south_north_dim
	1613	IF ( lati(i,j) .lt. 0. ) THEN
	1614	alpha(i,j) = - diff(i,j) * cone * RAD_PER_DEG
	1615	ELSE
	1616	alpha(i,j) = diff(i,j) * cone * RAD_PER_DEG
	1617	END IF
	1618	END DO
	1619	END DO
	1620
	1621
	1622	DO k = 1,bottom_top_dim
	1623	UUUmet(:,:,k) = VVV(:,:,k)sin(alpha) + UUU(:,:,k)cos(alpha)
	1624	VVVmet(:,:,k) = VVV(:,:,k)cos(alpha) - UUU(:,:,k)sin(alpha)
	1625	!print *,UUU(10,10,k), UUUmet(10,10,k)
	1626	!print *,VVV(10,10,k), VVVmet(10,10,k)
	1627	END DO
	1628	END IF
	1629
	1630	END SUBROUTINE calc_uvmet
	1631

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

Download in other formats:

Original Format