Context Navigation

gridcheck_nests.f90 @ 92fab65

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug

Last change on this file since 92fab65 was 92fab65, checked in by Ignacio Pisso <ip@…>, 5 years ago
add SPDX-License-Identifier to all .f90 files
Property mode set to `100644`
File size: 16.1 KB

Rev	Line
[92fab65]	1	! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
	2	! SPDX-License-Identifier: GPL-3.0-or-later
[332fbbd]	3
[e200b7a]	4	subroutine gridcheck_nests
	5
	6	!*****************************************************************************
	7	! *
	8	! This routine checks the grid specification for the nested model *
	9	! domains. It is similar to subroutine gridcheck, which checks the *
	10	! mother domain. *
	11	! *
	12	! Authors: A. Stohl, G. Wotawa *
	13	! *
	14	! 8 February 1999 *
	15	! *
	16	!*****************************************************************************
	17	! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api *
	18	! CHANGE: 03/12/2008, Harald Sodemann, change to f90 grib_api *
	19	!*****************************************************************************
	20
	21	use grib_api
	22	use par_mod
	23	use com_mod
	24
	25	implicit none
	26
	27	!HSO parameters for grib_api
	28	integer :: ifile
	29	integer :: iret
	30	integer :: igrib
	31	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
[4fbe7a5]	32	integer :: parID !added by mc for making it consistent with new gridcheck.f90
[e200b7a]	33	integer :: gotGrib
	34	!HSO end
	35	integer :: i,j,k,l,ifn,ifield,iumax,iwmax,numskip,nlev_ecn
	36	integer :: nuvzn,nwzn
	37	real :: akmn(nwzmax),bkmn(nwzmax),akzn(nuvzmax),bkzn(nuvzmax)
	38	real(kind=4) :: xaux1,xaux2,yaux1,yaux2
	39	real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
[4fbe7a5]	40	real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
[e200b7a]	41
	42	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
	43
	44	! dimension of isec2 at least (22+n), where n is the number of parallels or
	45	! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
	46
	47	! dimension of zsec2 at least (10+nn), where nn is the number of vertical
	48	! coordinate parameters
	49
	50	integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
	51	real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
	52
	53	!HSO grib api error messages
	54	character(len=24) :: gribErrorMsg = 'Error reading grib file'
	55	character(len=20) :: gribFunction = 'gridcheck_nests'
	56
	57	xresoln(0)=1. ! resolution enhancement for mother grid
	58	yresoln(0)=1. ! resolution enhancement for mother grid
	59
	60	! Loop about all nesting levels
	61	!******************************
	62
	63	do l=1,numbnests
	64
	65	iumax=0
	66	iwmax=0
	67
	68	if(ideltas.gt.0) then
	69	ifn=1
	70	else
	71	ifn=numbwf
	72	endif
	73	!
	74	! OPENING OF DATA FILE (GRIB CODE)
	75	!
	76	ifile=0
	77	igrib=0
	78	iret=0
	79
	80	5 call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
	81	(1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,ifn)),'r',iret)
	82	if (iret.ne.GRIB_SUCCESS) then
	83	goto 999 ! ERROR DETECTED
	84	endif
	85	!turn on support for multi fields messages
	86	!call grib_multi_support_on
	87
	88	gotGrib=0
	89	ifield=0
	90	10 ifield=ifield+1
	91
	92	!
	93	! GET NEXT FIELDS
	94	!
	95	call grib_new_from_file(ifile,igrib,iret)
	96	if (iret.eq.GRIB_END_OF_FILE) then
	97	goto 30 ! EOF DETECTED
	98	elseif (iret.ne.GRIB_SUCCESS) then
	99	goto 999 ! ERROR DETECTED
	100	endif
	101
	102	!first see if we read GRIB1 or GRIB2
	103	call grib_get_int(igrib,'editionNumber',gribVer,iret)
	104	call grib_check(iret,gribFunction,gribErrorMsg)
	105
	106	if (gribVer.eq.1) then ! GRIB Edition 1
	107
	108	!print*,'GRiB Edition 1'
	109	!read the grib2 identifiers
	110	call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
	111	call grib_check(iret,gribFunction,gribErrorMsg)
	112	call grib_get_int(igrib,'level',isec1(8),iret)
	113	call grib_check(iret,gribFunction,gribErrorMsg)
	114
	115	!change code for etadot to code for omega
	116	if (isec1(6).eq.77) then
	117	isec1(6)=135
	118	endif
	119
	120	!print*,isec1(6),isec1(8)
	121
	122	else
	123
	124	!print*,'GRiB Edition 2'
	125	!read the grib2 identifiers
	126	call grib_get_int(igrib,'discipline',discipl,iret)
	127	call grib_check(iret,gribFunction,gribErrorMsg)
	128	call grib_get_int(igrib,'parameterCategory',parCat,iret)
	129	call grib_check(iret,gribFunction,gribErrorMsg)
	130	call grib_get_int(igrib,'parameterNumber',parNum,iret)
	131	call grib_check(iret,gribFunction,gribErrorMsg)
	132	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
	133	call grib_check(iret,gribFunction,gribErrorMsg)
	134	call grib_get_int(igrib,'level',valSurf,iret)
	135	call grib_check(iret,gribFunction,gribErrorMsg)
[4fbe7a5]	136	call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new grid_check.f90
	137	call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new grid_check.f90
[e200b7a]	138
	139	!print*,discipl,parCat,parNum,typSurf,valSurf
	140
	141	!convert to grib1 identifiers
	142	isec1(6)=-1
	143	isec1(7)=-1
	144	isec1(8)=-1
	145	isec1(8)=valSurf ! level
	146	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
	147	isec1(6)=130 ! indicatorOfParameter
	148	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
	149	isec1(6)=131 ! indicatorOfParameter
	150	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
	151	isec1(6)=132 ! indicatorOfParameter
	152	elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
	153	isec1(6)=133 ! indicatorOfParameter
[b0434e1]	154	elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
	155	isec1(6)=246 ! indicatorOfParameter
	156	elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
	157	isec1(6)=247 ! indicatorOfParameter
	158	!ZHG end
	159	! ESO qc(=clwc+ciwc)
	160	elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
	161	isec1(6)=201031 ! indicatorOfParameter
[e200b7a]	162	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
	163	isec1(6)=134 ! indicatorOfParameter
	164	elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
	165	isec1(6)=135 ! indicatorOfParameter
[db712a8]	166	elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridcheck.f90
	167	isec1(6)=135 ! indicatorOfParameter !
[e200b7a]	168	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
	169	isec1(6)=151 ! indicatorOfParameter
	170	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
	171	isec1(6)=165 ! indicatorOfParameter
	172	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
	173	isec1(6)=166 ! indicatorOfParameter
	174	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
	175	isec1(6)=167 ! indicatorOfParameter
	176	elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
	177	isec1(6)=168 ! indicatorOfParameter
	178	elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
	179	isec1(6)=141 ! indicatorOfParameter
[4fbe7a5]	180	elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new gridchek.f90
[e200b7a]	181	isec1(6)=164 ! indicatorOfParameter
[4fbe7a5]	182	elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new gridchek.f90
[e200b7a]	183	isec1(6)=142 ! indicatorOfParameter
	184	elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
	185	isec1(6)=143 ! indicatorOfParameter
	186	elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
	187	isec1(6)=146 ! indicatorOfParameter
	188	elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
	189	isec1(6)=176 ! indicatorOfParameter
[4fbe7a5]	190	elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new gridchek.f90
[e200b7a]	191	isec1(6)=180 ! indicatorOfParameter
[4fbe7a5]	192	elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new gridchek.f90
[e200b7a]	193	isec1(6)=181 ! indicatorOfParameter
	194	elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
	195	isec1(6)=129 ! indicatorOfParameter
[4fbe7a5]	196	elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new gridchek.f90
[e200b7a]	197	isec1(6)=160 ! indicatorOfParameter
	198	elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
	199	(typSurf.eq.1)) then ! LSM
	200	isec1(6)=172 ! indicatorOfParameter
	201	else
	202	print,'**ERROR: undefined GRiB2 message found!',discipl, &
	203	parCat,parNum,typSurf
	204	endif
[4fbe7a5]	205	if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new gridchek.f90
	206	write(,) 'parId',parId, 'isec1(6)',isec1(6)
	207	! stop
	208	endif
[e200b7a]	209
	210	endif
	211
	212	!get the size and data of the values array
	213	if (isec1(6).ne.-1) then
	214	call grib_get_real4_array(igrib,'values',zsec4,iret)
	215	call grib_check(iret,gribFunction,gribErrorMsg)
	216	endif
	217
	218	!HSO get the required fields from section 2 in a gribex compatible manner
	219	if (ifield.eq.1) then
	220	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
	221	isec2(2),iret)
	222	call grib_check(iret,gribFunction,gribErrorMsg)
	223	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
	224	isec2(3),iret)
	225	call grib_check(iret,gribFunction,gribErrorMsg)
	226	call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
	227	isec2(12),iret)
	228	call grib_check(iret,gribFunction,gribErrorMsg)
	229	!HSO get the size and data of the vertical coordinate array
	230	call grib_get_real4_array(igrib,'pv',zsec2,iret)
	231	call grib_check(iret,gribFunction,gribErrorMsg)
	232
	233	nxn(l)=isec2(2)
	234	nyn(l)=isec2(3)
	235	nlev_ecn=isec2(12)/2-1
	236	endif ! ifield
	237
[db712a8]	238	if (nxn(l).gt.nxmaxn) then
	239	write(,) 'FLEXPART error: Too many grid points in x direction.'
	240	write(,) 'Reduce resolution of wind fields (file GRIDSPEC)'
	241	write(,) 'for nesting level ',l
	242	write(,) 'Or change parameter settings in file par_mod.'
	243	write(,) nxn(l),nxmaxn
	244	stop
	245	endif
	246
	247	if (nyn(l).gt.nymaxn) then
	248	write(,) 'FLEXPART error: Too many grid points in y direction.'
	249	write(,) 'Reduce resolution of wind fields (file GRIDSPEC)'
	250	write(,) 'for nesting level ',l
	251	write(,) 'Or change parameter settings in file par_mod.'
	252	write(,) nyn(l),nymaxn
	253	stop
	254	endif
	255
[e200b7a]	256	!HSO get the second part of the grid dimensions only from GRiB1 messages
[4fbe7a5]	257	if (isec1(6) .eq. 167 .and. (gotGrib.eq.0)) then !added by mc to make it consistent with new gridchek.f90 note that gotGrid must be changed in gotGrib!!
	258	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & !comment by mc: note that this was in the (if (ifield.eq.1) ..end above in gridchek.f90 see line 257
[e200b7a]	259	xaux1in,iret)
	260	call grib_check(iret,gribFunction,gribErrorMsg)
	261	call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
	262	xaux2in,iret)
	263	call grib_check(iret,gribFunction,gribErrorMsg)
	264	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
	265	yaux1in,iret)
	266	call grib_check(iret,gribFunction,gribErrorMsg)
	267	call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
	268	yaux2in,iret)
	269	call grib_check(iret,gribFunction,gribErrorMsg)
	270	xaux1=xaux1in
	271	xaux2=xaux2in
	272	yaux1=yaux1in
	273	yaux2=yaux2in
[4fbe7a5]	274	if(xaux1.gt.180.) xaux1=xaux1-360.0
	275	if(xaux2.gt.180.) xaux2=xaux2-360.0
	276	if(xaux1.lt.-180.) xaux1=xaux1+360.0
	277	if(xaux2.lt.-180.) xaux2=xaux2+360.0
	278	if (xaux2.lt.xaux1) xaux2=xaux2+360.0
[e200b7a]	279	xlon0n(l)=xaux1
	280	ylat0n(l)=yaux1
	281	dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
	282	dyn(l)=(yaux2-yaux1)/real(nyn(l)-1)
[4fbe7a5]	283	gotGrib=1 !commetn by mc note tahthere gotGRIB is used instead of gotGrid!!!
[e200b7a]	284	endif ! ifield.eq.1
	285
	286	k=isec1(8)
	287	if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
	288	if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
	289
	290	if(isec1(6).eq.129) then
	291	do j=0,nyn(l)-1
	292	do i=0,nxn(l)-1
	293	oron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
	294	end do
	295	end do
	296	endif
	297	if(isec1(6).eq.172) then
	298	do j=0,nyn(l)-1
	299	do i=0,nxn(l)-1
	300	lsmn(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
	301	end do
	302	end do
	303	endif
	304	if(isec1(6).eq.160) then
	305	do j=0,nyn(l)-1
	306	do i=0,nxn(l)-1
	307	excessoron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
	308	end do
	309	end do
	310	endif
	311
	312	call grib_release(igrib)
	313	goto 10 !! READ NEXT LEVEL OR PARAMETER
	314	!
	315	! CLOSING OF INPUT DATA FILE
	316	!
	317
	318	30 call grib_close_file(ifile)
	319
	320	!error message if no fields found with correct first longitude in it
	321	if (gotGrib.eq.0) then
	322	print,'**ERROR: input file needs to contain GRiB1 formatted'// &
	323	'messages'
	324	stop
	325	endif
	326
	327	nuvzn=iumax
	328	nwzn=iwmax
	329	if(nuvzn.eq.nlev_ec) nwzn=nlev_ecn+1
	330
	331	if ((nuvzn.gt.nuvzmax).or.(nwzn.gt.nwzmax)) then
	332	write(,) 'FLEXPART error: Nested wind fields have too many'// &
	333	'vertical levels.'
	334	write(,) 'Problem was encountered for nesting level ',l
	335	stop
	336	endif
	337
	338
	339	! Output of grid info
	340	!********************
	341
[b4d29ce]	342	write(*,'(a,i2,a)') ' Nested domain ',l,':'
[4fbe7a5]	343	write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Longitude range: ', &
[e200b7a]	344	xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
	345	' Grid distance: ',dxn(l)
[b4d29ce]	346	write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Latitude range : ', &
[e200b7a]	347	ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
	348	' Grid distance: ',dyn(l)
	349	write(,)
	350
	351	! Determine, how much the resolutions in the nests are enhanced as
	352	! compared to the mother grid
	353	!*****************************************************************
	354
	355	xresoln(l)=dx/dxn(l)
	356	yresoln(l)=dy/dyn(l)
	357
	358	! Determine the mother grid coordinates of the corner points of the
	359	! nested grids
	360	! Convert first to geographical coordinates, then to grid coordinates
	361	!********************************************************************
	362
	363	xaux1=xlon0n(l)
	364	xaux2=xlon0n(l)+real(nxn(l)-1)*dxn(l)
	365	yaux1=ylat0n(l)
	366	yaux2=ylat0n(l)+real(nyn(l)-1)*dyn(l)
	367
	368	xln(l)=(xaux1-xlon0)/dx
	369	xrn(l)=(xaux2-xlon0)/dx
	370	yln(l)=(yaux1-ylat0)/dy
	371	yrn(l)=(yaux2-ylat0)/dy
	372
	373
	374	if ((xln(l).lt.0.).or.(yln(l).lt.0.).or. &
	375	(xrn(l).gt.real(nxmin1)).or.(yrn(l).gt.real(nymin1))) then
	376	write(,) 'Nested domain does not fit into mother domain'
	377	write(,) 'For global mother domain fields, you can shift'
	378	write(,) 'shift the mother domain into x-direction'
	379	write(,) 'by setting nxshift (file par_mod) to a'
	380	write(,) 'positive value. Execution is terminated.'
	381	stop
	382	endif
	383
	384
	385	! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
	386	! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
	387
	388	numskip=nlev_ecn-nuvzn ! number of ecmwf model layers not used by FLEXPART
	389	do i=1,nwzn
	390	j=numskip+i
	391	k=nlev_ecn+1+numskip+i
	392	akmn(nwzn-i+1)=zsec2(j)
	393	bkmn(nwzn-i+1)=zsec2(k)
	394	end do
	395
	396	!
	397	! CALCULATION OF AKZ, BKZ
	398	! AKZ,BKZ: model discretization parameters at the center of each model
	399	! layer
	400	!
	401	! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
	402	! i.e. ground level
	403	!*****************************************************************************
	404
	405	akzn(1)=0.
	406	bkzn(1)=1.0
	407	do i=1,nuvzn
	408	akzn(i+1)=0.5*(akmn(i+1)+akmn(i))
	409	bkzn(i+1)=0.5*(bkmn(i+1)+bkmn(i))
	410	end do
	411	nuvzn=nuvzn+1
	412
	413	! Check, whether the heights of the model levels of the nested
	414	! wind fields are consistent with those of the mother domain.
	415	! If not, terminate model run.
	416	!*************************************************************
	417
	418	do i=1,nuvz
	419	if ((akzn(i).ne.akz(i)).or.(bkzn(i).ne.bkz(i))) then
	420	write(,) 'FLEXPART error: The wind fields of nesting level',l
	421	write(,) 'are not consistent with the mother domain:'
	422	write(,) 'Differences in vertical levels detected.'
	423	stop
	424	endif
	425	end do
	426
	427	do i=1,nwz
	428	if ((akmn(i).ne.akm(i)).or.(bkmn(i).ne.bkm(i))) then
	429	write(,) 'FLEXPART error: The wind fields of nesting level',l
	430	write(,) 'are not consistent with the mother domain:'
	431	write(,) 'Differences in vertical levels detected.'
	432	stop
	433	endif
	434	end do
	435
	436	end do
	437
	438	return
	439
	440	999 write(,)
	441	write(,) ' ###########################################'// &
	442	'###### '
	443	write(,) ' FLEXPART SUBROUTINE GRIDCHECK:'
	444	write(,) ' CAN NOT OPEN INPUT DATA FILE '//wfnamen(l,ifn)
	445	write(,) ' FOR NESTING LEVEL ',k
	446	write(,) ' ###########################################'// &
	447	'###### '
	448	stop
	449
	450	end subroutine gridcheck_nests

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

Download in other formats:

Original Format