Context Navigation

gridcheck.f90 @ 4

Last change on this file since 4 was 4, checked in by mlanger, 11 years ago

Property svn:executable set to ``*
File size: 20.2 KB

Rev	Line
[4]	1	!**********************************************************************
	2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
	3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
	4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
	5	! *
	6	! This file is part of FLEXPART. *
	7	! *
	8	! FLEXPART is free software: you can redistribute it and/or modify *
	9	! it under the terms of the GNU General Public License as published by*
	10	! the Free Software Foundation, either version 3 of the License, or *
	11	! (at your option) any later version. *
	12	! *
	13	! FLEXPART is distributed in the hope that it will be useful, *
	14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
	15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	16	! GNU General Public License for more details. *
	17	! *
	18	! You should have received a copy of the GNU General Public License *
	19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
	20	!**********************************************************************
	21
	22	subroutine gridcheck
	23
	24	!**********************************************************************
	25	! *
	26	! FLEXPART MODEL SUBROUTINE GRIDCHECK *
	27	! *
	28	!**********************************************************************
	29	! *
	30	! AUTHOR: G. WOTAWA *
	31	! DATE: 1997-08-06 *
	32	! LAST UPDATE: 1997-10-10 *
	33	! *
	34	! Update: 1999-02-08, global fields allowed, A. Stohl*
	35	! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
	36	! ECMWF grib_api *
	37	! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
	38	! ECMWF grib_api *
	39	! *
	40	!**********************************************************************
	41	! *
	42	! DESCRIPTION: *
	43	! *
	44	! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT *
	45	! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST- *
	46	! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE *
	47	! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES *
	48	! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT *
	49	! ANY CALL. *
	50	! *
	51	! XLON0 geographical longitude of lower left gridpoint *
	52	! YLAT0 geographical latitude of lower left gridpoint *
	53	! NX number of grid points x-direction *
	54	! NY number of grid points y-direction *
	55	! DX grid distance x-direction *
	56	! DY grid distance y-direction *
	57	! NUVZ number of grid points for horizontal wind *
	58	! components in z direction *
	59	! NWZ number of grid points for vertical wind *
	60	! component in z direction *
	61	! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
	62	! points of the polar stereographic grid): *
	63	! used to check the CFL criterion *
	64	! UVHEIGHT(1)- heights of gridpoints where u and v are *
	65	! UVHEIGHT(NUVZ) given *
	66	! WHEIGHT(1)- heights of gridpoints where w is given *
	67	! WHEIGHT(NWZ) *
	68	! *
	69	!**********************************************************************
	70
	71	use grib_api
	72	use par_mod
	73	use com_mod
	74	use conv_mod
	75	use cmapf_mod, only: stlmbr,stcm2p
	76
	77	implicit none
	78
	79	!HSO parameters for grib_api
	80	integer :: ifile
	81	integer :: iret
	82	integer :: igrib
	83	integer :: gotGrid
	84	real(kind=4) :: xaux1,xaux2,yaux1,yaux2
	85	real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
	86	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
	87	!HSO end
	88	integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
	89	real :: sizesouth,sizenorth,xauxa,pint
	90
	91	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
	92
	93	! dimension of isec2 at least (22+n), where n is the number of parallels or
	94	! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
	95
	96	! dimension of zsec2 at least (10+nn), where nn is the number of vertical
	97	! coordinate parameters
	98
	99	integer :: isec1(56),isec2(22+nxmax+nymax)
	100	real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
	101	character(len=1) :: opt
	102
	103	!HSO grib api error messages
	104	character(len=24) :: gribErrorMsg = 'Error reading grib file'
	105	character(len=20) :: gribFunction = 'gridcheck'
	106
	107	iumax=0
	108	iwmax=0
	109
	110	if(ideltas.gt.0) then
	111	ifn=1
	112	else
	113	ifn=numbwf
	114	endif
	115	!
	116	! OPENING OF DATA FILE (GRIB CODE)
	117	!
	118	5 call grib_open_file(ifile,path(3)(1:length(3)) &
	119	//trim(wfname(ifn)),'r',iret)
	120	if (iret.ne.GRIB_SUCCESS) then
	121	goto 999 ! ERROR DETECTED
	122	endif
	123	!turn on support for multi fields messages
	124	!call grib_multi_support_on
	125
	126	gotGrid=0
	127	ifield=0
	128	10 ifield=ifield+1
	129
	130	!
	131	! GET NEXT FIELDS
	132	!
	133	call grib_new_from_file(ifile,igrib,iret)
	134	if (iret.eq.GRIB_END_OF_FILE ) then
	135	goto 30 ! EOF DETECTED
	136	elseif (iret.ne.GRIB_SUCCESS) then
	137	goto 999 ! ERROR DETECTED
	138	endif
	139
	140	!first see if we read GRIB1 or GRIB2
	141	call grib_get_int(igrib,'editionNumber',gribVer,iret)
	142	call grib_check(iret,gribFunction,gribErrorMsg)
	143
	144	if (gribVer.eq.1) then ! GRIB Edition 1
	145
	146	!print*,'GRiB Edition 1'
	147	!read the grib2 identifiers
	148	call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
	149	call grib_check(iret,gribFunction,gribErrorMsg)
	150	call grib_get_int(igrib,'level',isec1(8),iret)
	151	call grib_check(iret,gribFunction,gribErrorMsg)
	152
	153	!change code for etadot to code for omega
	154	if (isec1(6).eq.77) then
	155	isec1(6)=135
	156	endif
	157
	158	!print*,isec1(6),isec1(8)
	159
	160	else
	161
	162	!print*,'GRiB Edition 2'
	163	!read the grib2 identifiers
	164	call grib_get_int(igrib,'discipline',discipl,iret)
	165	call grib_check(iret,gribFunction,gribErrorMsg)
	166	call grib_get_int(igrib,'parameterCategory',parCat,iret)
	167	call grib_check(iret,gribFunction,gribErrorMsg)
	168	call grib_get_int(igrib,'parameterNumber',parNum,iret)
	169	call grib_check(iret,gribFunction,gribErrorMsg)
	170	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
	171	call grib_check(iret,gribFunction,gribErrorMsg)
	172	call grib_get_int(igrib,'level',valSurf,iret)
	173	call grib_check(iret,gribFunction,gribErrorMsg)
	174
	175	!print*,discipl,parCat,parNum,typSurf,valSurf
	176
	177	!convert to grib1 identifiers
	178	isec1(6)=-1
	179	isec1(7)=-1
	180	isec1(8)=-1
	181	isec1(8)=valSurf ! level
	182	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
	183	isec1(6)=130 ! indicatorOfParameter
	184	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
	185	isec1(6)=131 ! indicatorOfParameter
	186	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
	187	isec1(6)=132 ! indicatorOfParameter
	188	elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
	189	isec1(6)=133 ! indicatorOfParameter
	190	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
	191	isec1(6)=134 ! indicatorOfParameter
	192	elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
	193	isec1(6)=135 ! indicatorOfParameter
	194	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
	195	isec1(6)=151 ! indicatorOfParameter
	196	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
	197	isec1(6)=165 ! indicatorOfParameter
	198	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
	199	isec1(6)=166 ! indicatorOfParameter
	200	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
	201	isec1(6)=167 ! indicatorOfParameter
	202	elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
	203	isec1(6)=168 ! indicatorOfParameter
	204	elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
	205	isec1(6)=141 ! indicatorOfParameter
	206	elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
	207	isec1(6)=164 ! indicatorOfParameter
	208	elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
	209	isec1(6)=142 ! indicatorOfParameter
	210	elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
	211	isec1(6)=143 ! indicatorOfParameter
	212	elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
	213	isec1(6)=146 ! indicatorOfParameter
	214	elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
	215	isec1(6)=176 ! indicatorOfParameter
	216	elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
	217	isec1(6)=180 ! indicatorOfParameter
	218	elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
	219	isec1(6)=181 ! indicatorOfParameter
	220	elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
	221	isec1(6)=129 ! indicatorOfParameter
	222	elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
	223	isec1(6)=160 ! indicatorOfParameter
	224	elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
	225	(typSurf.eq.1)) then ! LSM
	226	isec1(6)=172 ! indicatorOfParameter
	227	else
	228	print,'**ERROR: undefined GRiB2 message found!',discipl, &
	229	parCat,parNum,typSurf
	230	endif
	231
	232	endif
	233
	234	!get the size and data of the values array
	235	if (isec1(6).ne.-1) then
	236	call grib_get_real4_array(igrib,'values',zsec4,iret)
	237	call grib_check(iret,gribFunction,gribErrorMsg)
	238	endif
	239
	240	if (ifield.eq.1) then
	241
	242	!HSO get the required fields from section 2 in a gribex compatible manner
	243	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
	244	isec2(2),iret)
	245	call grib_check(iret,gribFunction,gribErrorMsg)
	246	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
	247	isec2(3),iret)
	248	call grib_check(iret,gribFunction,gribErrorMsg)
	249	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
	250	xaux1in,iret)
	251	call grib_check(iret,gribFunction,gribErrorMsg)
	252	call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
	253	isec2(12),iret)
	254	call grib_check(iret,gribFunction,gribErrorMsg)
	255
	256	! get the size and data of the vertical coordinate array
	257	call grib_get_real4_array(igrib,'pv',zsec2,iret)
	258	call grib_check(iret,gribFunction,gribErrorMsg)
	259
	260	nxfield=isec2(2)
	261	ny=isec2(3)
	262	nlev_ec=isec2(12)/2-1
	263	endif
	264
	265	!HSO get the second part of the grid dimensions only from GRiB1 messages
	266	if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
	267	call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
	268	xaux2in,iret)
	269	call grib_check(iret,gribFunction,gribErrorMsg)
	270	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
	271	yaux1in,iret)
	272	call grib_check(iret,gribFunction,gribErrorMsg)
	273	call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
	274	yaux2in,iret)
	275	call grib_check(iret,gribFunction,gribErrorMsg)
	276	xaux1=xaux1in
	277	xaux2=xaux2in
	278	yaux1=yaux1in
	279	yaux2=yaux2in
	280	if (xaux1.gt.180.) xaux1=xaux1-360.0
	281	if (xaux2.gt.180.) xaux2=xaux2-360.0
	282	if (xaux1.lt.-180.) xaux1=xaux1+360.0
	283	if (xaux2.lt.-180.) xaux2=xaux2+360.0
	284	if (xaux2.lt.xaux1) xaux2=xaux2+360.0
	285	xlon0=xaux1
	286	ylat0=yaux1
	287	dx=(xaux2-xaux1)/real(nxfield-1)
	288	dy=(yaux2-yaux1)/real(ny-1)
	289	dxconst=180./(dxr_earthpi)
	290	dyconst=180./(dyr_earthpi)
	291	gotGrid=1
	292
	293	! Check whether fields are global
	294	! If they contain the poles, specify polar stereographic map
	295	! projections using the stlmbr- and stcm2p-calls
	296	!***********************************************************
	297
	298	xauxa=abs(xaux2+dx-360.-xaux1)
	299	if (xauxa.lt.0.001) then
	300	nx=nxfield+1 ! field is cyclic
	301	xglobal=.true.
	302	if (abs(nxshift).ge.nx) &
	303	stop 'nxshift in file par_mod is too large'
	304	xlon0=xlon0+real(nxshift)*dx
	305	else
	306	nx=nxfield
	307	xglobal=.false.
	308	if (nxshift.ne.0) &
	309	stop 'nxshift (par_mod) must be zero for non-global domain'
	310	endif
	311	nxmin1=nx-1
	312	nymin1=ny-1
	313	if (xlon0.gt.180.) xlon0=xlon0-360.
	314	xauxa=abs(yaux1+90.)
	315	if (xglobal.and.xauxa.lt.0.001) then
	316	sglobal=.true. ! field contains south pole
	317	! Enhance the map scale by factor 3 (*2=6) compared to north-south
	318	! map scale
	319	sizesouth=6.*(switchsouth+90.)/dy
	320	call stlmbr(southpolemap,-90.,0.)
	321	call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
	322	sizesouth,switchsouth,180.)
	323	switchsouthg=(switchsouth-ylat0)/dy
	324	else
	325	sglobal=.false.
	326	switchsouthg=999999.
	327	endif
	328	xauxa=abs(yaux2-90.)
	329	if (xglobal.and.xauxa.lt.0.001) then
	330	nglobal=.true. ! field contains north pole
	331	! Enhance the map scale by factor 3 (*2=6) compared to north-south
	332	! map scale
	333	sizenorth=6.*(90.-switchnorth)/dy
	334	call stlmbr(northpolemap,90.,0.)
	335	call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
	336	sizenorth,switchnorth,180.)
	337	switchnorthg=(switchnorth-ylat0)/dy
	338	else
	339	nglobal=.false.
	340	switchnorthg=999999.
	341	endif
	342	if (nxshift.lt.0) &
	343	stop 'nxshift (par_mod) must not be negative'
	344	if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
	345	endif ! gotGrid
	346
	347	k=isec1(8)
	348	if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
	349	if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
	350
	351	if(isec1(6).eq.129) then
	352	do jy=0,ny-1
	353	do ix=0,nxfield-1
	354	oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
	355	end do
	356	end do
	357	endif
	358	if(isec1(6).eq.172) then
	359	do jy=0,ny-1
	360	do ix=0,nxfield-1
	361	lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
	362	end do
	363	end do
	364	endif
	365	if(isec1(6).eq.160) then
	366	do jy=0,ny-1
	367	do ix=0,nxfield-1
	368	excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
	369	end do
	370	end do
	371	endif
	372
	373	call grib_release(igrib)
	374	goto 10 !! READ NEXT LEVEL OR PARAMETER
	375	!
	376	! CLOSING OF INPUT DATA FILE
	377	!
	378
	379	30 call grib_close_file(ifile)
	380
	381	!error message if no fields found with correct first longitude in it
	382	if (gotGrid.eq.0) then
	383	print,'**ERROR: input file needs to contain GRiB1 formatted'// &
	384	'messages'
	385	stop
	386	endif
	387
	388	nuvz=iumax
	389	nwz =iwmax
	390	if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
	391
	392	if (nx.gt.nxmax) then
	393	write(,) 'FLEXPART error: Too many grid points in x direction.'
	394	write(,) 'Reduce resolution of wind fields.'
	395	write(,) 'Or change parameter settings in file par_mod.'
	396	write(,) nx,nxmax
	397	stop
	398	endif
	399
	400	if (ny.gt.nymax) then
	401	write(,) 'FLEXPART error: Too many grid points in y direction.'
	402	write(,) 'Reduce resolution of wind fields.'
	403	write(,) 'Or change parameter settings in file par_mod.'
	404	write(,) ny,nymax
	405	stop
	406	endif
	407
	408	if (nuvz+1.gt.nuvzmax) then
	409	write(,) 'FLEXPART error: Too many u,v grid points in z '// &
	410	'direction.'
	411	write(,) 'Reduce resolution of wind fields.'
	412	write(,) 'Or change parameter settings in file par_mod.'
	413	write(,) nuvz+1,nuvzmax
	414	stop
	415	endif
	416
	417	if (nwz.gt.nwzmax) then
	418	write(,) 'FLEXPART error: Too many w grid points in z '// &
	419	'direction.'
	420	write(,) 'Reduce resolution of wind fields.'
	421	write(,) 'Or change parameter settings in file par_mod.'
	422	write(,) nwz,nwzmax
	423	stop
	424	endif
	425
	426	! If desired, shift all grids by nxshift grid cells
	427	!**************************************************
	428
	429	if (xglobal) then
	430	call shift_field_0(oro,nxfield,ny)
	431	call shift_field_0(lsm,nxfield,ny)
	432	call shift_field_0(excessoro,nxfield,ny)
	433	endif
	434
	435	! Output of grid info
	436	!********************
	437
	438	write(,)
	439	write(,)
	440	write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
	441	nuvz+1,nwz
	442	write(,)
	443	write(*,'(a)') 'Mother domain:'
	444	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
	445	xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
	446	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ', &
	447	ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
	448	write(,)
	449
	450
	451	! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
	452	! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
	453
	454	numskip=nlev_ec-nuvz ! number of ecmwf model layers not used
	455	! by trajectory model
	456	!do 8940 i=1,244
	457	! write (,) 'zsec2:',i,ifield,zsec2(i),numskip
	458	!940 continue
	459	! stop
	460	! SEC SEC SEC
	461	! for unknown reason zsec 1 to 10 is filled in this version
	462	! compared to the old one
	463	! SEC SEC SE
	464	do i=1,nwz
	465	j=numskip+i
	466	k=nlev_ec+1+numskip+i
	467	akm(nwz-i+1)=zsec2(j)
	468	! write (,) 'ifield:',ifield,k,j,zsec2(10+j)
	469	bkm(nwz-i+1)=zsec2(k)
	470	end do
	471
	472	!
	473	! CALCULATION OF AKZ, BKZ
	474	! AKZ,BKZ: model discretization parameters at the center of each model
	475	! layer
	476	!
	477	! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
	478	! i.e. ground level
	479	!*****************************************************************************
	480
	481	akz(1)=0.
	482	bkz(1)=1.0
	483	do i=1,nuvz
	484	akz(i+1)=0.5*(akm(i+1)+akm(i))
	485	bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
	486	end do
	487	nuvz=nuvz+1
	488
	489	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
	490	! upon the transformation to z levels. In order to save computer memory, this is
	491	! not done anymore in the standard version. However, this option can still be
	492	! switched on by replacing the following lines with those below, that are
	493	! currently commented out. For this, similar changes are necessary in
	494	! verttransform.f and verttranform_nests.f
	495	!*****************************************************************************
	496
	497	nz=nuvz
	498	if (nz.gt.nzmax) stop 'nzmax too small'
	499	do i=1,nuvz
	500	aknew(i)=akz(i)
	501	bknew(i)=bkz(i)
	502	end do
	503
	504	! Switch on following lines to use doubled vertical resolution
	505	!*************************************************************
	506	!nz=nuvz+nwz-1
	507	!if (nz.gt.nzmax) stop 'nzmax too small'
	508	!do 100 i=1,nwz
	509	! aknew(2*(i-1)+1)=akm(i)
	510	!00 bknew(2*(i-1)+1)=bkm(i)
	511	!do 110 i=2,nuvz
	512	! aknew(2*(i-1))=akz(i)
	513	!10 bknew(2*(i-1))=bkz(i)
	514	! End doubled vertical resolution
	515
	516
	517	! Determine the uppermost level for which the convection scheme shall be applied
	518	! by assuming that there is no convection above 50 hPa (for standard SLP)
	519	!*****************************************************************************
	520
	521	do i=1,nuvz-2
	522	pint=akz(i)+bkz(i)*101325.
	523	if (pint.lt.5000.) goto 96
	524	end do
	525	96 nconvlev=i
	526	if (nconvlev.gt.nconvlevmax-1) then
	527	nconvlev=nconvlevmax-1
	528	write(,) 'Attention, convection only calculated up to ', &
	529	akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
	530	endif
	531
	532	return
	533
	534	999 write(,)
	535	write(,) ' ###########################################'// &
	536	'###### '
	537	write(,) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
	538	write(,) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
	539	write(,) ' ###########################################'// &
	540	'###### '
	541	write(,)
	542	write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND !!!'
	543	write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE... !!!'
	544	write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
	545	write(*,'(a)') '!!! THE "X" KEY... !!!'
	546	write(,)
	547	read(*,'(a)') opt
	548	if(opt.eq.'X') then
	549	stop
	550	else
	551	goto 5
	552	endif
	553
	554	end subroutine gridcheck

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

Download in other formats:

Original Format