Context Navigation

readwind_gfs.f90 @ 9ca6e38

GFS_025dev

Last change on this file since 9ca6e38 was 9ca6e38, checked in by Espen Sollum ATMOS <eso@…>, 3 years ago
minor edits
Property mode set to `100755`
File size: 25.1 KB

Rev	Line
[92fab65]	1	! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
	2	! SPDX-License-Identifier: GPL-3.0-or-later
[e200b7a]	3
[6ecb30a]	4	subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
[e200b7a]	5
[9ca6e38]	6	!***********************************************************************
	7	!* *
	8	!* TRAJECTORY MODEL SUBROUTINE READWIND *
	9	!* *
	10	!***********************************************************************
	11	!* *
	12	!* AUTHOR: G. WOTAWA *
	13	!* DATE: 1997-08-05 *
	14	!* LAST UPDATE: 2000-10-17, Andreas Stohl *
	15	!* CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
	16	!* qvh (on eta coordinates) in common block *
	17	!* CHANGE: 16/11/2005, Caroline Forster, GFS data *
	18	!* CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2 *
	19	!* data with the ECMWF grib_api library *
	20	!* CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
	21	!* ECMWF grib_api *
	22	! *
	23	! Unified ECMWF and GFS builds *
	24	! Marian Harustak, 12.5.2017 *
	25	! - Renamed routine from readwind to readwind_gfs *
	26	!* *
	27	!***********************************************************************
	28	!* *
	29	!* DESCRIPTION: *
	30	!* *
	31	!* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE *
	32	!* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE *
	33	!* *
	34	!* INPUT: *
	35	!* indj indicates number of the wind field to be read in *
	36	!* n temporal index for meteorological fields (1 to 3)*
	37	!* *
	38	!* IMPORTANT VARIABLES FROM COMMON BLOCK: *
	39	!* *
	40	!* wfname File name of data to be read in *
	41	!* nx,ny,nuvz,nwz expected field dimensions *
	42	!* nlev_ec number of vertical levels ecmwf model *
	43	!* uu,vv,ww wind fields *
	44	!* tt,qv temperature and specific humidity *
	45	!* ps surface pressure *
	46	!* *
	47	!***********************************************************************
[e200b7a]	48
[a756649]	49	use eccodes
[e200b7a]	50	use par_mod
	51	use com_mod
	52
	53	implicit none
	54
[9ca6e38]	55	!HSO new parameters for grib_api
[e200b7a]	56	integer :: ifile
	57	integer :: iret
	58	integer :: igrib
[9ca6e38]	59	integer :: gribVer,parCat,parNum,typSurf,discipl,valSurf
	60	!HSO end edits
[e200b7a]	61	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
	62	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
	63	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
	64	integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
	65
[9ca6e38]	66	! NCEP
[5d74ed9]	67	integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
[e200b7a]	68	real :: help, temp, ew
	69	real :: elev
	70	real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
	71	real :: tlev1(0:nxmax-1,0:nymax-1)
	72	real :: qvh2(0:nxmax-1,0:nymax-1)
	73
[467460a]	74	integer :: i180
[e200b7a]	75
[9ca6e38]	76	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
	77	!HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
[e200b7a]	78
	79	integer :: isec1(8),isec2(3)
[467460a]	80	real :: xsec18
[e200b7a]	81	real(kind=4) :: zsec4(jpunp)
	82	real(kind=4) :: xaux,yaux,xaux0,yaux0
[467460a]	83	real,parameter :: eps=spacing(2.0_4*360.0_4)
[e200b7a]	84	real(kind=8) :: xauxin,yauxin
	85	real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
	86	real :: plev1,hlev1,ff10m,fflev1
	87
	88	logical :: hflswitch,strswitch
	89
[9ca6e38]	90	!HSO for grib api error messages
[e200b7a]	91	character(len=24) :: gribErrorMsg = 'Error reading grib file'
	92	character(len=20) :: gribFunction = 'readwind_gfs'
[db91eb7]	93	character(len=20) :: shortname
[e200b7a]	94
	95
	96	hflswitch=.false.
	97	strswitch=.false.
	98	levdiff2=nlev_ec-nwz+1
	99	iumax=0
	100	iwmax=0
	101
	102
[9ca6e38]	103	! OPENING OF DATA FILE (GRIB CODE)
[e200b7a]	104
[9ca6e38]	105	!HSO
[db91eb7]	106	call grib_open_file(ifile,path(3)(1:length(3)) &
[9ca6e38]	107	//trim(wfname(indj)),'r',iret)
[e200b7a]	108	if (iret.ne.GRIB_SUCCESS) then
	109	goto 888 ! ERROR DETECTED
	110	endif
[9ca6e38]	111	!turn on support for multi fields messages
[e200b7a]	112	call grib_multi_support_on
	113
	114	numpt=0
	115	numpu=0
	116	numpv=0
	117	numpw=0
	118	numprh=0
[7873bf7]	119	numpclwch=0
[e200b7a]	120	ifield=0
[9ca6e38]	121	10 ifield=ifield+1
	122	!
	123	! GET NEXT FIELDS
	124	!
[e200b7a]	125	call grib_new_from_file(ifile,igrib,iret)
	126	if (iret.eq.GRIB_END_OF_FILE) then
	127	goto 50 ! EOF DETECTED
	128	elseif (iret.ne.GRIB_SUCCESS) then
	129	goto 888 ! ERROR DETECTED
	130	endif
	131
[9ca6e38]	132	!first see if we read GRIB1 or GRIB2
[e200b7a]	133	call grib_get_int(igrib,'editionNumber',gribVer,iret)
[8a65cb0]	134	! call grib_check(iret,gribFunction,gribErrorMsg)
[e200b7a]	135
	136	if (gribVer.eq.1) then ! GRIB Edition 1
	137
[9ca6e38]	138	!read the grib1 identifiers
	139
	140	call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
	141	call grib_check(iret,gribFunction,gribErrorMsg)
	142	call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
	143	call grib_check(iret,gribFunction,gribErrorMsg)
	144	call grib_get_int(igrib,'level',isec1(8),iret)
	145	call grib_check(iret,gribFunction,gribErrorMsg)
	146	!JMA / SH: isec1(8) not evaluated any more below
	147	!b/c with GRIB 2 this may be a real variable
	148	xsec18 = real(isec1(8))
[e200b7a]	149
	150	else ! GRIB Edition 2
	151
[9ca6e38]	152	!read the grib2 identifiers
	153
	154	call grib_get_string(igrib,'shortName',shortname,iret)
	155
	156	call grib_get_int(igrib,'discipline',discipl,iret)
	157	call grib_check(iret,gribFunction,gribErrorMsg)
	158	call grib_get_int(igrib,'parameterCategory',parCat,iret)
	159	call grib_check(iret,gribFunction,gribErrorMsg)
	160	call grib_get_int(igrib,'parameterNumber',parNum,iret)
	161	call grib_check(iret,gribFunction,gribErrorMsg)
	162	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
	163	call grib_check(iret,gribFunction,gribErrorMsg)
	164	!
	165	call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
	166	valSurf,iret)
	167	call grib_check(iret,gribFunction,gribErrorMsg)
	168
	169	! write(,) 'Field: ',ifield,parCat,parNum,typSurf,shortname
	170	!convert to grib1 identifiers
	171	isec1(6:8)=-1
	172	xsec18 =-1.0
	173	!JMA / SH: isec1(8) not evaluated any more below
	174	!b/c with GRIB 2 this may be a real variable
	175	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
	176	isec1(6)=11 ! indicatorOfParameter
	177	isec1(7)=100 ! indicatorOfTypeOfLevel
	178	xsec18=valSurf/100.0 ! level, convert to hPa
	179	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
	180	isec1(6)=33 ! indicatorOfParameter
	181	isec1(7)=100 ! indicatorOfTypeOfLevel
	182	xsec18=valSurf/100.0 ! level, convert to hPa
	183	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
	184	isec1(6)=34 ! indicatorOfParameter
	185	isec1(7)=100 ! indicatorOfTypeOfLevel
	186	xsec18=valSurf/100.0 ! level, convert to hPa
	187	elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
	188	isec1(6)=39 ! indicatorOfParameter
	189	isec1(7)=100 ! indicatorOfTypeOfLevel
	190	xsec18=valSurf/100.0 ! level, convert to hPa
	191	elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
	192	isec1(6)=52 ! indicatorOfParameter
	193	isec1(7)=100 ! indicatorOfTypeOfLevel
	194	xsec18=valSurf/100.0 ! level, convert to hPa
	195	elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
	196	isec1(6)=52 ! indicatorOfParameter
	197	isec1(7)=105 ! indicatorOfTypeOfLevel
	198	xsec18=real(2)
	199	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
	200	isec1(6)=11 ! indicatorOfParameter
	201	isec1(7)=105 ! indicatorOfTypeOfLevel
	202	xsec18=real(2)
	203	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
	204	isec1(6)=33 ! indicatorOfParameter
	205	isec1(7)=105 ! indicatorOfTypeOfLevel
	206	xsec18=real(10)
	207	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
	208	isec1(6)=34 ! indicatorOfParameter
	209	isec1(7)=105 ! indicatorOfTypeOfLevel
	210	xsec18=real(10)
	211	elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
	212	isec1(6)=153 ! indicatorOfParameter
	213	isec1(7)=100 ! indicatorOfTypeOfLevel
	214	xsec18=valSurf/100.0 ! level, convert to hPa
	215	elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
	216	isec1(6)=2 ! indicatorOfParameter
	217	isec1(7)=102 ! indicatorOfTypeOfLevel
	218	xsec18=real(0)
	219	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
	220	isec1(6)=1 ! indicatorOfParameter
	221	isec1(7)=1 ! indicatorOfTypeOfLevel
	222	xsec18=real(0)
	223	elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
	224	isec1(6)=66 ! indicatorOfParameter
	225	isec1(7)=1 ! indicatorOfTypeOfLevel
	226	xsec18=real(0)
	227	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
	228	isec1(6)=11 ! indicatorOfParameter
	229	isec1(7)=107 ! indicatorOfTypeOfLevel
	230	xsec18=0.995 ! lowest sigma level
	231	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
	232	isec1(6)=33 ! indicatorOfParameter
	233	isec1(7)=107 ! indicatorOfTypeOfLevel
	234	xsec18=0.995 ! lowest sigma level
	235	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
	236	isec1(6)=34 ! indicatorOfParameter
	237	isec1(7)=107 ! indicatorOfTypeOfLevel
	238	xsec18=0.995 ! lowest sigma level
	239	elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
	240	isec1(6)=7 ! indicatorOfParameter
	241	isec1(7)=1 ! indicatorOfTypeOfLevel
	242	xsec18=real(0)
	243	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
	244	.and.(discipl.eq.2)) then ! LSM
	245	isec1(6)=81 ! indicatorOfParameter
	246	isec1(7)=1 ! indicatorOfTypeOfLevel
	247	xsec18=real(0)
	248	elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
	249	isec1(6)=221 ! indicatorOfParameter
	250	isec1(7)=1 ! indicatorOfTypeOfLevel
	251	xsec18=real(0)
	252	elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
	253	isec1(6)=62 ! indicatorOfParameter
	254	isec1(7)=1 ! indicatorOfTypeOfLevel
	255	xsec18=real(0)
	256	elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
	257	isec1(6)=63 ! indicatorOfParameter
	258	isec1(7)=1 ! indicatorOfTypeOfLevel
	259	xsec18=real(0)
	260	endif
[e200b7a]	261
	262	endif ! gribVer
	263
	264	if (isec1(6).ne.-1) then
[9ca6e38]	265	! get the size and data of the values array
[e200b7a]	266	call grib_get_real4_array(igrib,'values',zsec4,iret)
[467460a]	267	call grib_check(iret,gribFunction,gribErrorMsg)
[e200b7a]	268	endif
	269
	270	if(ifield.eq.1) then
	271
[9ca6e38]	272	!get the required fields from section 2
	273	!store compatible to gribex input
	274	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
	275	isec2(2),iret)
	276	call grib_check(iret,gribFunction,gribErrorMsg)
	277	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
	278	isec2(3),iret)
	279	call grib_check(iret,gribFunction,gribErrorMsg)
	280	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
	281	xauxin,iret)
	282	call grib_check(iret,gribFunction,gribErrorMsg)
	283	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
	284	yauxin,iret)
	285	call grib_check(iret,gribFunction,gribErrorMsg)
	286	xaux=xauxin+real(nxshift)*dx
	287	yaux=yauxin
	288
	289	! CHECK GRID SPECIFICATIONS
[e200b7a]	290
	291	if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
	292	if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
[467460a]	293	if(xaux.eq.0.) xaux=-180.0 ! NCEP DATA
[e200b7a]	294	xaux0=xlon0
	295	yaux0=ylat0
	296	if(xaux.lt.0.) xaux=xaux+360.
	297	if(yaux.lt.0.) yaux=yaux+360.
	298	if(xaux0.lt.0.) xaux0=xaux0+360.
	299	if(yaux0.lt.0.) yaux0=yaux0+360.
[467460a]	300	if(abs(xaux-xaux0).gt.eps) then
[9ca6e38]	301	write (, ) xaux, xaux0
	302	stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
[467460a]	303	endif
	304	if(abs(yaux-yaux0).gt.eps) then
[9ca6e38]	305	write (, ) yaux, yaux0
	306	stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
[467460a]	307	end if
[e200b7a]	308	endif
[9ca6e38]	309	!HSO end of edits
[e200b7a]	310
[467460a]	311	i180=nint(180./dx)
[e200b7a]	312
	313	if (isec1(6).ne.-1) then
	314
[467460a]	315	! write (, ) 'nxfield: ', nxfield, i180
	316
[9ca6e38]	317	do j=0,nymin1
	318	do i=0,nxfield-1
	319	if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
[467460a]	320	! TEMPERATURE
[9ca6e38]	321	if((i.eq.0).and.(j.eq.0)) then
	322	numpt=minloc(abs(xsec18*100.0-akz),dim=1)
	323	endif
	324	help=zsec4(nxfield*(ny-j-1)+i+1)
	325	if (help.eq.0) then
	326	write (, ) 'i, j: ', i, j
	327	stop 'help == 0.0'
	328	endif
	329	if(i.lt.i180) then
	330	tth(i180+i,j,numpt,n)=help
	331	else
	332	tth(i-i180,j,numpt,n)=help
	333	endif
[e200b7a]	334	endif
[9ca6e38]	335	if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
	336	! U VELOCITY
	337	if((i.eq.0).and.(j.eq.0)) then
	338	numpu=minloc(abs(xsec18*100.0-akz),dim=1)
	339	endif
	340	help=zsec4(nxfield*(ny-j-1)+i+1)
	341	if(i.lt.i180) then
	342	uuh(i180+i,j,numpu)=help
	343	else
	344	uuh(i-i180,j,numpu)=help
	345	endif
[e200b7a]	346	endif
[9ca6e38]	347	if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
	348	! V VELOCITY
	349	if((i.eq.0).and.(j.eq.0)) then
	350	numpv=minloc(abs(xsec18*100.0-akz),dim=1)
	351	endif
	352	help=zsec4(nxfield*(ny-j-1)+i+1)
	353	if(i.lt.i180) then
	354	vvh(i180+i,j,numpv)=help
	355	else
	356	vvh(i-i180,j,numpv)=help
	357	endif
[e200b7a]	358	endif
[9ca6e38]	359	if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
	360	! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
	361	if((i.eq.0).and.(j.eq.0)) then
	362	numprh=minloc(abs(xsec18*100.0-akz),dim=1)
	363	endif
	364	help=zsec4(nxfield*(ny-j-1)+i+1)
	365	if(i.lt.i180) then
	366	qvh(i180+i,j,numprh,n)=help
	367	else
	368	qvh(i-i180,j,numprh,n)=help
	369	endif
[e200b7a]	370	endif
[9ca6e38]	371	if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
	372	! SURFACE PRESSURE
	373	help=zsec4(nxfield*(ny-j-1)+i+1)
	374	if(i.lt.i180) then
	375	ps(i180+i,j,1,n)=help
	376	else
	377	ps(i-i180,j,1,n)=help
	378	endif
[e200b7a]	379	endif
[9ca6e38]	380	if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
	381	! W VELOCITY
	382	if((i.eq.0).and.(j.eq.0)) numpw=minloc(abs(xsec18*100.0-akz),dim=1)
	383	help=zsec4(nxfield*(ny-j-1)+i+1)
	384	if(i.lt.i180) then
	385	wwh(i180+i,j,numpw)=help
	386	else
	387	wwh(i-i180,j,numpw)=help
	388	endif
[e200b7a]	389	endif
[9ca6e38]	390	if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
	391	! SNOW DEPTH
	392	help=zsec4(nxfield*(ny-j-1)+i+1)
	393	if(i.lt.i180) then
	394	sd(i180+i,j,1,n)=help
	395	else
	396	sd(i-i180,j,1,n)=help
	397	endif
[e200b7a]	398	endif
[9ca6e38]	399	if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
	400	! MEAN SEA LEVEL PRESSURE
	401	help=zsec4(nxfield*(ny-j-1)+i+1)
	402	if(i.lt.i180) then
	403	msl(i180+i,j,1,n)=help
	404	else
	405	msl(i-i180,j,1,n)=help
	406	endif
[e200b7a]	407	endif
[9ca6e38]	408	if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
	409	! TOTAL CLOUD COVER
	410	help=zsec4(nxfield*(ny-j-1)+i+1)
	411	if(i.lt.i180) then
	412	tcc(i180+i,j,1,n)=help
	413	else
	414	tcc(i-i180,j,1,n)=help
	415	endif
[e200b7a]	416	endif
[9ca6e38]	417	if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
	418	(nint(xsec18).eq.10)) then
	419	! 10 M U VELOCITY
	420	help=zsec4(nxfield*(ny-j-1)+i+1)
	421	if(i.lt.i180) then
	422	u10(i180+i,j,1,n)=help
	423	else
	424	u10(i-i180,j,1,n)=help
	425	endif
[e200b7a]	426	endif
[9ca6e38]	427	if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
	428	(nint(xsec18).eq.10)) then
	429	! 10 M V VELOCITY
	430	help=zsec4(nxfield*(ny-j-1)+i+1)
	431	if(i.lt.i180) then
	432	v10(i180+i,j,1,n)=help
	433	else
	434	v10(i-i180,j,1,n)=help
	435	endif
[e200b7a]	436	endif
[9ca6e38]	437	if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
	438	(nint(xsec18).eq.2)) then
	439	! 2 M TEMPERATURE
	440	help=zsec4(nxfield*(ny-j-1)+i+1)
	441	if(i.lt.i180) then
	442	tt2(i180+i,j,1,n)=help
	443	else
	444	tt2(i-i180,j,1,n)=help
	445	endif
[e200b7a]	446	endif
[9ca6e38]	447	if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
	448	(nint(xsec18).eq.2)) then
	449	! 2 M DEW POINT TEMPERATURE
	450	help=zsec4(nxfield*(ny-j-1)+i+1)
	451	if(i.lt.i180) then
	452	td2(i180+i,j,1,n)=help
	453	else
	454	td2(i-i180,j,1,n)=help
	455	endif
[e200b7a]	456	endif
[9ca6e38]	457	if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
	458	! LARGE SCALE PREC.
	459	help=zsec4(nxfield*(ny-j-1)+i+1)
	460	if(i.lt.i180) then
	461	lsprec(i180+i,j,1,n)=help
	462	else
	463	lsprec(i-i180,j,1,n)=help
	464	endif
[e200b7a]	465	endif
[9ca6e38]	466	if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
	467	! CONVECTIVE PREC.
	468	help=zsec4(nxfield*(ny-j-1)+i+1)
	469	if(i.lt.i180) then
	470	convprec(i180+i,j,1,n)=help
	471	else
	472	convprec(i-i180,j,1,n)=help
	473	endif
[e200b7a]	474	endif
[9ca6e38]	475	if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
	476	! TOPOGRAPHY
	477	help=zsec4(nxfield*(ny-j-1)+i+1)
	478	if(i.lt.i180) then
	479	oro(i180+i,j)=help
	480	excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
	481	else
	482	oro(i-i180,j)=help
	483	excessoro(i-i180,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
	484	endif
[e200b7a]	485	endif
[9ca6e38]	486	if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
	487	! LAND SEA MASK
	488	help=zsec4(nxfield*(ny-j-1)+i+1)
	489	if(i.lt.i180) then
	490	lsm(i180+i,j)=help
	491	else
	492	lsm(i-i180,j)=help
	493	endif
[e200b7a]	494	endif
[9ca6e38]	495	if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
	496	! MIXING HEIGHT
	497	help=zsec4(nxfield*(ny-j-1)+i+1)
	498	if(i.lt.i180) then
	499	hmix(i180+i,j,1,n)=help
	500	else
	501	hmix(i-i180,j,1,n)=help
	502	endif
[e200b7a]	503	endif
[9ca6e38]	504	if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
	505	(nint(xsec18).eq.02)) then
	506	! 2 M RELATIVE HUMIDITY
	507	help=zsec4(nxfield*(ny-j-1)+i+1)
	508	if(i.lt.i180) then
	509	qvh2(i180+i,j)=help
	510	else
	511	qvh2(i-i180,j)=help
	512	endif
[e200b7a]	513	endif
[9ca6e38]	514	if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
	515	! TEMPERATURE LOWEST SIGMA LEVEL
	516	help=zsec4(nxfield*(ny-j-1)+i+1)
	517	if(i.lt.i180) then
	518	tlev1(i180+i,j)=help
	519	else
	520	tlev1(i-i180,j)=help
	521	endif
[e200b7a]	522	endif
[9ca6e38]	523	if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
	524	! U VELOCITY LOWEST SIGMA LEVEL
	525	help=zsec4(nxfield*(ny-j-1)+i+1)
	526	if(i.lt.i180) then
	527	ulev1(i180+i,j)=help
	528	else
	529	ulev1(i-i180,j)=help
	530	endif
[e200b7a]	531	endif
[9ca6e38]	532	if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
	533	! V VELOCITY LOWEST SIGMA LEVEL
	534	help=zsec4(nxfield*(ny-j-1)+i+1)
	535	if(i.lt.i180) then
	536	vlev1(i180+i,j)=help
	537	else
	538	vlev1(i-i180,j)=help
	539	endif
[e200b7a]	540	endif
[db91eb7]	541	! SEC & IP 12/2018 read GFS clouds
[9ca6e38]	542	if((isec1(6).eq.153).and.(isec1(7).eq.100)) then !! CLWCR Cloud liquid water content [kg/kg]
	543	if((i.eq.0).and.(j.eq.0)) then
	544	numpclwch=minloc(abs(xsec18*100.0-akz),dim=1)
	545	endif
	546	help=zsec4(nxfield*(ny-j-1)+i+1)
	547	if(i.lt.i180) then
	548	clwch(i180+i,j,numpclwch,n)=help
	549	else
	550	clwch(i-i180,j,numpclwch,n)=help
	551	endif
	552	readclouds=.true.
	553	sumclouds=.true.
	554	endif
[db91eb7]	555
[e200b7a]	556
[9ca6e38]	557	end do
[e200b7a]	558	end do
	559
	560	endif
	561
	562	if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
[9ca6e38]	563	! NCEP ISOBARIC LEVELS
[e200b7a]	564	iumax=iumax+1
	565	endif
	566
	567	call grib_release(igrib)
	568	goto 10 !! READ NEXT LEVEL OR PARAMETER
[9ca6e38]	569	!
	570	! CLOSING OF INPUT DATA FILE
	571	!
[e200b7a]	572
[9ca6e38]	573	!HSO close grib file
	574	50 continue
[e200b7a]	575	call grib_close_file(ifile)
	576
[9ca6e38]	577	! SENS. HEAT FLUX
[e200b7a]	578	sshf(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files
	579	hflswitch=.false. ! Heat flux not available
[9ca6e38]	580	! SOLAR RADIATIVE FLUXES
[e200b7a]	581	ssr(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files
[9ca6e38]	582	! EW SURFACE STRESS
[e200b7a]	583	ewss=0.0 ! not available from gfs.tccz.pgrbfxx files
[9ca6e38]	584	! NS SURFACE STRESS
[e200b7a]	585	nsss=0.0 ! not available from gfs.tccz.pgrbfxx files
	586	strswitch=.false. ! stress not available
	587
[9ca6e38]	588	! CONVERT TP TO LSP (GRIB2 only)
[e200b7a]	589	if (gribVer.eq.2) then
[467460a]	590	lsprec(0:nxfield-1, 0:nymin1, 1, n) = max( 0.0 , lsprec(0:nxfield-1,0:nymin1,1,n)-convprec(0:nxfield-1,0:nymin1,1,n) )
[e200b7a]	591	endif
[9ca6e38]	592	!HSO end edits
[e200b7a]	593
	594
[9ca6e38]	595	! TRANSFORM RH TO SPECIFIC HUMIDITY
[e200b7a]	596
	597	do j=0,ny-1
	598	do i=0,nxfield-1
	599	do k=1,nuvz
	600	help=qvh(i,j,k,n)
	601	temp=tth(i,j,k,n)
[467460a]	602	if (temp .eq. 0.0) then
[9ca6e38]	603	write (, ) i, j, k, n
	604	temp = 273.0
[467460a]	605	endif
[e200b7a]	606	plev1=akm(k)+bkm(k)*ps(i,j,1,n)
	607	elev=ew(temp)*help/100.0
	608	qvh(i,j,k,n)=xmwml(elev/(plev1-((1.0-xmwml)elev)))
	609	end do
	610	end do
	611	end do
	612
[9ca6e38]	613	! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
	614	! USING BOLTON'S (1980) FORMULA
	615	! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
[e200b7a]	616
	617	do j=0,ny-1
	618	do i=0,nxfield-1
[9ca6e38]	619	help=qvh2(i,j)
	620	temp=tt2(i,j,1,n)
	621	if (temp .eq. 0.0) then
	622	write (, ) i, j, n
	623	temp = 273.0
	624	endif
	625	elev=ew(temp)/100.*help/100. !vapour pressure in hPa
	626	td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
	627	if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
[e200b7a]	628	end do
	629	end do
	630
	631	if(levdiff2.eq.0) then
	632	iwmax=nlev_ec+1
	633	do i=0,nxmin1
	634	do j=0,nymin1
	635	wwh(i,j,nlev_ec+1)=0.
	636	end do
	637	end do
	638	endif
	639
	640
[9ca6e38]	641	! For global fields, assign the leftmost data column also to the rightmost
	642	! data column; if required, shift whole grid by nxshift grid points
	643	!*************************************************************************
[e200b7a]	644
	645	if (xglobal) then
	646	call shift_field_0(ewss,nxfield,ny)
	647	call shift_field_0(nsss,nxfield,ny)
	648	call shift_field_0(oro,nxfield,ny)
	649	call shift_field_0(excessoro,nxfield,ny)
	650	call shift_field_0(lsm,nxfield,ny)
	651	call shift_field_0(ulev1,nxfield,ny)
	652	call shift_field_0(vlev1,nxfield,ny)
	653	call shift_field_0(tlev1,nxfield,ny)
	654	call shift_field_0(qvh2,nxfield,ny)
	655	call shift_field(ps,nxfield,ny,1,1,2,n)
	656	call shift_field(sd,nxfield,ny,1,1,2,n)
	657	call shift_field(msl,nxfield,ny,1,1,2,n)
	658	call shift_field(tcc,nxfield,ny,1,1,2,n)
	659	call shift_field(u10,nxfield,ny,1,1,2,n)
	660	call shift_field(v10,nxfield,ny,1,1,2,n)
	661	call shift_field(tt2,nxfield,ny,1,1,2,n)
	662	call shift_field(td2,nxfield,ny,1,1,2,n)
	663	call shift_field(lsprec,nxfield,ny,1,1,2,n)
	664	call shift_field(convprec,nxfield,ny,1,1,2,n)
	665	call shift_field(sshf,nxfield,ny,1,1,2,n)
	666	call shift_field(ssr,nxfield,ny,1,1,2,n)
	667	call shift_field(hmix,nxfield,ny,1,1,2,n)
	668	call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
	669	call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
	670	call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
	671	call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
	672	call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
[db91eb7]	673	! IP & SEC adding GFS Clouds 20181205
	674	call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
[e200b7a]	675	endif
	676
	677	do i=0,nxmin1
	678	do j=0,nymin1
[9ca6e38]	679	! Convert precip. from mm/s -> mm/hour
[e200b7a]	680	convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
	681	lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
	682	surfstr(i,j,1,n)=sqrt(ewss(i,j)2+nsss(i,j)2)
	683	end do
	684	end do
	685
	686	if ((.not.hflswitch).or.(.not.strswitch)) then
[9ca6e38]	687	! write(,) 'WARNING: No flux data contained in GRIB file ',
	688	! + wfname(indj)
[e200b7a]	689
[9ca6e38]	690	! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
	691	!***************************************************************************
[e200b7a]	692
	693	do i=0,nxmin1
	694	do j=0,nymin1
	695	hlev1=30.0 ! HEIGHT OF FIRST MODEL SIGMA LAYER
	696	ff10m= sqrt(u10(i,j,1,n)2+v10(i,j,1,n)2)
	697	fflev1=sqrt(ulev1(i,j)2+vlev1(i,j)2)
	698	call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
	699	tt2(i,j,1,n),tlev1(i,j),ff10m,fflev1, &
	700	surfstr(i,j,1,n),sshf(i,j,1,n))
	701	if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
	702	if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
	703	end do
	704	end do
	705	endif
	706
	707	if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
	708	if(iumax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT'
	709
	710	return
[9ca6e38]	711	888 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
[e200b7a]	712	write(,) ' #### ',wfname(indj),' #### '
	713	write(,) ' #### IS NOT GRIB FORMAT !!! #### '
	714	stop 'Execution terminated'
[9ca6e38]	715	999 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
[e200b7a]	716	write(,) ' #### ',wfname(indj),' #### '
	717	write(,) ' #### CANNOT BE OPENED !!! #### '
	718	stop 'Execution terminated'
	719
[6ecb30a]	720	end subroutine readwind_gfs

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

Download in other formats:

Original Format