Context Navigation

gridcheck.f90 @ 20

Last change on this file since 20 was 20, checked in by igpis, 10 years ago
move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others
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)
[20]	100	!real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
	101	real(kind=4) :: zsec2(184),zsec4(jpunp)
[4]	102	character(len=1) :: opt
	103
	104	!HSO grib api error messages
	105	character(len=24) :: gribErrorMsg = 'Error reading grib file'
	106	character(len=20) :: gribFunction = 'gridcheck'
	107
	108	iumax=0
	109	iwmax=0
	110
	111	if(ideltas.gt.0) then
	112	ifn=1
	113	else
	114	ifn=numbwf
	115	endif
	116	!
	117	! OPENING OF DATA FILE (GRIB CODE)
	118	!
	119	5 call grib_open_file(ifile,path(3)(1:length(3)) &
	120	//trim(wfname(ifn)),'r',iret)
	121	if (iret.ne.GRIB_SUCCESS) then
	122	goto 999 ! ERROR DETECTED
	123	endif
	124	!turn on support for multi fields messages
	125	!call grib_multi_support_on
	126
	127	gotGrid=0
	128	ifield=0
	129	10 ifield=ifield+1
	130
	131	!
	132	! GET NEXT FIELDS
	133	!
	134	call grib_new_from_file(ifile,igrib,iret)
	135	if (iret.eq.GRIB_END_OF_FILE ) then
	136	goto 30 ! EOF DETECTED
	137	elseif (iret.ne.GRIB_SUCCESS) then
	138	goto 999 ! ERROR DETECTED
	139	endif
	140
	141	!first see if we read GRIB1 or GRIB2
	142	call grib_get_int(igrib,'editionNumber',gribVer,iret)
	143	call grib_check(iret,gribFunction,gribErrorMsg)
	144
	145	if (gribVer.eq.1) then ! GRIB Edition 1
	146
	147	!print*,'GRiB Edition 1'
	148	!read the grib2 identifiers
	149	call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
	150	call grib_check(iret,gribFunction,gribErrorMsg)
	151	call grib_get_int(igrib,'level',isec1(8),iret)
	152	call grib_check(iret,gribFunction,gribErrorMsg)
	153
	154	!change code for etadot to code for omega
	155	if (isec1(6).eq.77) then
	156	isec1(6)=135
	157	endif
	158
	159	!print*,isec1(6),isec1(8)
	160
	161	else
	162
	163	!print*,'GRiB Edition 2'
	164	!read the grib2 identifiers
	165	call grib_get_int(igrib,'discipline',discipl,iret)
	166	call grib_check(iret,gribFunction,gribErrorMsg)
	167	call grib_get_int(igrib,'parameterCategory',parCat,iret)
	168	call grib_check(iret,gribFunction,gribErrorMsg)
	169	call grib_get_int(igrib,'parameterNumber',parNum,iret)
	170	call grib_check(iret,gribFunction,gribErrorMsg)
	171	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
	172	call grib_check(iret,gribFunction,gribErrorMsg)
	173	call grib_get_int(igrib,'level',valSurf,iret)
	174	call grib_check(iret,gribFunction,gribErrorMsg)
	175
	176	!print*,discipl,parCat,parNum,typSurf,valSurf
	177
	178	!convert to grib1 identifiers
	179	isec1(6)=-1
	180	isec1(7)=-1
	181	isec1(8)=-1
	182	isec1(8)=valSurf ! level
	183	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
	184	isec1(6)=130 ! indicatorOfParameter
	185	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
	186	isec1(6)=131 ! indicatorOfParameter
	187	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
	188	isec1(6)=132 ! indicatorOfParameter
	189	elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
	190	isec1(6)=133 ! indicatorOfParameter
	191	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
	192	isec1(6)=134 ! indicatorOfParameter
	193	elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
	194	isec1(6)=135 ! indicatorOfParameter
	195	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
	196	isec1(6)=151 ! indicatorOfParameter
	197	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
	198	isec1(6)=165 ! indicatorOfParameter
	199	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
	200	isec1(6)=166 ! indicatorOfParameter
	201	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
	202	isec1(6)=167 ! indicatorOfParameter
	203	elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
	204	isec1(6)=168 ! indicatorOfParameter
	205	elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
	206	isec1(6)=141 ! indicatorOfParameter
	207	elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
	208	isec1(6)=164 ! indicatorOfParameter
	209	elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
	210	isec1(6)=142 ! indicatorOfParameter
	211	elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
	212	isec1(6)=143 ! indicatorOfParameter
	213	elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
	214	isec1(6)=146 ! indicatorOfParameter
	215	elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
	216	isec1(6)=176 ! indicatorOfParameter
	217	elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
	218	isec1(6)=180 ! indicatorOfParameter
	219	elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
	220	isec1(6)=181 ! indicatorOfParameter
	221	elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
	222	isec1(6)=129 ! indicatorOfParameter
	223	elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
	224	isec1(6)=160 ! indicatorOfParameter
	225	elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
	226	(typSurf.eq.1)) then ! LSM
	227	isec1(6)=172 ! indicatorOfParameter
	228	else
	229	print,'**ERROR: undefined GRiB2 message found!',discipl, &
	230	parCat,parNum,typSurf
	231	endif
	232
	233	endif
	234
	235	!get the size and data of the values array
	236	if (isec1(6).ne.-1) then
	237	call grib_get_real4_array(igrib,'values',zsec4,iret)
	238	call grib_check(iret,gribFunction,gribErrorMsg)
	239	endif
	240
	241	if (ifield.eq.1) then
	242
	243	!HSO get the required fields from section 2 in a gribex compatible manner
	244	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
	245	isec2(2),iret)
	246	call grib_check(iret,gribFunction,gribErrorMsg)
	247	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
	248	isec2(3),iret)
	249	call grib_check(iret,gribFunction,gribErrorMsg)
	250	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
	251	xaux1in,iret)
	252	call grib_check(iret,gribFunction,gribErrorMsg)
	253	call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
	254	isec2(12),iret)
	255	call grib_check(iret,gribFunction,gribErrorMsg)
	256
	257	! get the size and data of the vertical coordinate array
	258	call grib_get_real4_array(igrib,'pv',zsec2,iret)
	259	call grib_check(iret,gribFunction,gribErrorMsg)
	260
	261	nxfield=isec2(2)
	262	ny=isec2(3)
	263	nlev_ec=isec2(12)/2-1
	264	endif
	265
	266	!HSO get the second part of the grid dimensions only from GRiB1 messages
	267	if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
	268	call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
	269	xaux2in,iret)
	270	call grib_check(iret,gribFunction,gribErrorMsg)
	271	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
	272	yaux1in,iret)
	273	call grib_check(iret,gribFunction,gribErrorMsg)
	274	call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
	275	yaux2in,iret)
	276	call grib_check(iret,gribFunction,gribErrorMsg)
	277	xaux1=xaux1in
	278	xaux2=xaux2in
	279	yaux1=yaux1in
	280	yaux2=yaux2in
	281	if (xaux1.gt.180.) xaux1=xaux1-360.0
	282	if (xaux2.gt.180.) xaux2=xaux2-360.0
	283	if (xaux1.lt.-180.) xaux1=xaux1+360.0
	284	if (xaux2.lt.-180.) xaux2=xaux2+360.0
	285	if (xaux2.lt.xaux1) xaux2=xaux2+360.0
	286	xlon0=xaux1
	287	ylat0=yaux1
	288	dx=(xaux2-xaux1)/real(nxfield-1)
	289	dy=(yaux2-yaux1)/real(ny-1)
	290	dxconst=180./(dxr_earthpi)
	291	dyconst=180./(dyr_earthpi)
	292	gotGrid=1
	293
	294	! Check whether fields are global
	295	! If they contain the poles, specify polar stereographic map
	296	! projections using the stlmbr- and stcm2p-calls
	297	!***********************************************************
	298
	299	xauxa=abs(xaux2+dx-360.-xaux1)
	300	if (xauxa.lt.0.001) then
	301	nx=nxfield+1 ! field is cyclic
	302	xglobal=.true.
	303	if (abs(nxshift).ge.nx) &
	304	stop 'nxshift in file par_mod is too large'
	305	xlon0=xlon0+real(nxshift)*dx
	306	else
	307	nx=nxfield
	308	xglobal=.false.
	309	if (nxshift.ne.0) &
	310	stop 'nxshift (par_mod) must be zero for non-global domain'
	311	endif
	312	nxmin1=nx-1
	313	nymin1=ny-1
	314	if (xlon0.gt.180.) xlon0=xlon0-360.
	315	xauxa=abs(yaux1+90.)
	316	if (xglobal.and.xauxa.lt.0.001) then
	317	sglobal=.true. ! field contains south pole
	318	! Enhance the map scale by factor 3 (*2=6) compared to north-south
	319	! map scale
	320	sizesouth=6.*(switchsouth+90.)/dy
	321	call stlmbr(southpolemap,-90.,0.)
	322	call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
	323	sizesouth,switchsouth,180.)
	324	switchsouthg=(switchsouth-ylat0)/dy
	325	else
	326	sglobal=.false.
	327	switchsouthg=999999.
	328	endif
	329	xauxa=abs(yaux2-90.)
	330	if (xglobal.and.xauxa.lt.0.001) then
	331	nglobal=.true. ! field contains north pole
	332	! Enhance the map scale by factor 3 (*2=6) compared to north-south
	333	! map scale
	334	sizenorth=6.*(90.-switchnorth)/dy
	335	call stlmbr(northpolemap,90.,0.)
	336	call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
	337	sizenorth,switchnorth,180.)
	338	switchnorthg=(switchnorth-ylat0)/dy
	339	else
	340	nglobal=.false.
	341	switchnorthg=999999.
	342	endif
	343	if (nxshift.lt.0) &
	344	stop 'nxshift (par_mod) must not be negative'
	345	if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
	346	endif ! gotGrid
	347
	348	k=isec1(8)
	349	if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
	350	if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
	351
	352	if(isec1(6).eq.129) then
	353	do jy=0,ny-1
	354	do ix=0,nxfield-1
	355	oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
	356	end do
	357	end do
	358	endif
	359	if(isec1(6).eq.172) then
	360	do jy=0,ny-1
	361	do ix=0,nxfield-1
	362	lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
	363	end do
	364	end do
	365	endif
	366	if(isec1(6).eq.160) then
	367	do jy=0,ny-1
	368	do ix=0,nxfield-1
	369	excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
	370	end do
	371	end do
	372	endif
	373
	374	call grib_release(igrib)
	375	goto 10 !! READ NEXT LEVEL OR PARAMETER
	376	!
	377	! CLOSING OF INPUT DATA FILE
	378	!
	379
	380	30 call grib_close_file(ifile)
	381
	382	!error message if no fields found with correct first longitude in it
	383	if (gotGrid.eq.0) then
	384	print,'**ERROR: input file needs to contain GRiB1 formatted'// &
	385	'messages'
	386	stop
	387	endif
	388
	389	nuvz=iumax
	390	nwz =iwmax
	391	if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
	392
	393	if (nx.gt.nxmax) then
	394	write(,) 'FLEXPART error: Too many grid points in x direction.'
	395	write(,) 'Reduce resolution of wind fields.'
	396	write(,) 'Or change parameter settings in file par_mod.'
	397	write(,) nx,nxmax
	398	stop
	399	endif
	400
	401	if (ny.gt.nymax) then
	402	write(,) 'FLEXPART error: Too many grid points in y direction.'
	403	write(,) 'Reduce resolution of wind fields.'
	404	write(,) 'Or change parameter settings in file par_mod.'
	405	write(,) ny,nymax
	406	stop
	407	endif
	408
	409	if (nuvz+1.gt.nuvzmax) then
	410	write(,) 'FLEXPART error: Too many u,v grid points in z '// &
	411	'direction.'
	412	write(,) 'Reduce resolution of wind fields.'
	413	write(,) 'Or change parameter settings in file par_mod.'
	414	write(,) nuvz+1,nuvzmax
	415	stop
	416	endif
	417
	418	if (nwz.gt.nwzmax) then
	419	write(,) 'FLEXPART error: Too many w grid points in z '// &
	420	'direction.'
	421	write(,) 'Reduce resolution of wind fields.'
	422	write(,) 'Or change parameter settings in file par_mod.'
	423	write(,) nwz,nwzmax
	424	stop
	425	endif
	426
	427	! If desired, shift all grids by nxshift grid cells
	428	!**************************************************
	429
	430	if (xglobal) then
	431	call shift_field_0(oro,nxfield,ny)
	432	call shift_field_0(lsm,nxfield,ny)
	433	call shift_field_0(excessoro,nxfield,ny)
	434	endif
	435
	436	! Output of grid info
	437	!********************
	438
	439	write(,)
	440	write(,)
	441	write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
	442	nuvz+1,nwz
	443	write(,)
	444	write(*,'(a)') 'Mother domain:'
	445	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
	446	xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
	447	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ', &
	448	ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
	449	write(,)
	450
	451
	452	! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
	453	! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
	454
	455	numskip=nlev_ec-nuvz ! number of ecmwf model layers not used
	456	! by trajectory model
	457	!do 8940 i=1,244
	458	! write (,) 'zsec2:',i,ifield,zsec2(i),numskip
	459	!940 continue
	460	! stop
	461	! SEC SEC SEC
	462	! for unknown reason zsec 1 to 10 is filled in this version
	463	! compared to the old one
	464	! SEC SEC SE
	465	do i=1,nwz
	466	j=numskip+i
	467	k=nlev_ec+1+numskip+i
	468	akm(nwz-i+1)=zsec2(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