Context Navigation

readwind.f90 @ 5f2e8f6

flexpart-noresm

Last change on this file since 5f2e8f6 was 5f2e8f6, checked in by Ignacio Pisso <Ignacio.Pisso@…>, 8 years ago
new flexpart noresm code in src
Property mode set to `100755`
File size: 31.2 KB

Line
1	!**********************************************************************
2	! Copyright 2016 *
3	! Andreas Stohl, Massimo Cassiani, Petra Seibert, A. Frank, *
4	! Gerhard Wotawa, Caroline Forster, Sabine Eckhardt, John Burkhart, *
5	! Harald Sodemann, Ignacio Pisso *
6	! *
7	! This file is part of FLEXPART-NorESM *
8	! *
9	! FLEXPART-NorESM is free software: you can redistribute it *
10	! and/or modify *
11	! it under the terms of the GNU General Public License as published by*
12	! the Free Software Foundation, either version 3 of the License, or *
13	! (at your option) any later version. *
14	! *
15	! FLEXPART-NorESM is distributed in the hope that it will be useful, *
16	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	! GNU General Public License for more details. *
19	! *
20	! You should have received a copy of the GNU General Public License *
21	! along with FLEXPART-NorESM. *
22	! If not, see <http://www.gnu.org/licenses/>. *
23	!**********************************************************************
24
25
26
27	Subroutine readwind(indj,n,uuh,vvh,wwh)
28
29
30	!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
31	!* *
32	!* Read all fields needed to run FLEXPART *
33	!* note: part of netcdf reading strucutre adapted from the routines by *
34	!* Fast and Easter in FLEXPART-WRF *
35	!* if methodw is 1 it call a ruotine to obtain etadotdpdeta *
36	!* *
37	!* *
38	!* Author: *
39	!* M. Cassiani 2016 *
40	!* *
41	!* *
42	!* *
43	!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
44
45
46
47	use noresm_variables
48	use par_mod
49	use com_mod
50	use conv_mod
51	use cmapf_mod, only: stlmbr,stcm2p
52
53	implicit none
54
55	include 'netcdf.inc'
56
57	!integer maxdim,maxvar,maxtime !moved to par_mod.f90
58	real(kind=4) :: eps
59
60	!parameter(maxdim=9,maxvar=60,maxtime=12, !moved in in common block par_mod
61	parameter(eps=0.0001)
62
63	integer :: gotGrid
64	real(kind=4) :: sizesouth,sizenorth,xauxa,pint
65
66	!c---------------------------------------------
67	real(kind=4) :: xaux1,xaux2,yaux1,yaux2
68	real(kind=dp) :: xaux1in,xaux2in,yaux1in,yaux2in
69
70	integer :: nyfieldin,nxfieldin
71
72	!C-------------- dichiarazioni per netcdf
73	real(kind=4) :: dewpoint ! this is a function
74	real(kind=4) :: duma,twomdewpoint
75	real(kind=4), allocatable, dimension(:) :: duma_alloc
76	integer :: date_aid(maxtime),idumb,itime ,dimtimenum,indextime, &
77	time_interval
78	real(kind=dp) :: jul,juldate,dumb !variable jul used for date in days, juldate is a function
79
80	integer :: ierr !error code message
81	integer :: idiagaa !flag
82	integer :: id_var,id_dim(maxdim)
83	integer :: nvar_exp_in_grid_atm_nc, &
84	nvar_exp_in_meteo_field
85
86
87	integer :: i,ii,j,k, iatt, idimid_unlim, idum, iret, ivtype
88	integer :: ix,jy,kz
89	integer :: lenatt, lendim(maxdim)
90	integer :: natts_tot, ncid, ndims_tot, nvars_tot
91	integer :: sizetype
92	character*110 :: fnamenc
93	character*80 :: dimname(maxdim)
94	character*80 :: attname
95	character*160 :: varname,vartype
96	character*160 :: varnamev(maxvar)
97	character*160 :: units(maxvar)
98	character*160 :: vartypev(maxvar)
99	integer :: ndimsv(maxvar)
100	integer :: dimidsv(maxvar,maxdim)
101	character*1000 :: dumch1000
102
103	!integer istart(maxdim),icount(maxdim)
104	integer :: xtype,xtypev(maxvar)
105	integer :: ndims
106	integer :: dimids(maxdim),dimidsm(maxvar,maxdim)
107	integer :: LENDIM_EXP(maxdim),LENDIM_MAX(maxdim)
108	integer :: natts
109	integer :: varid
110
111	real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
112	real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
113	real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
114	real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),&
115	ewss(0:nxmax-1,0:nymax-1)
116	real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
117
118	logical :: hflswitch,strswitch
119	integer :: indj,n !i,j,k,n,levdiff2,ifield,iumax,iwmax
120	character*7 :: stringwftime
121
122
123	!***************************************************************************************************
124	9100 format( / '** read_noresmout_gridinfo -- ', a /&
125	'file = ', a )
126	9110 format( / '** read_noresmout_gridinfo -- ', a, 1x, i8 /&
127	'file = ', a )
128	9120 format( / '** read_noresmout_gridinfo -- ', a, 2(1x,i8) /&
129	'file = ', a )
130
131	9030 format( a, 2i6, 2(2x,a) )
132	9031 format( 1i4,1a20,2i4,1a30)
133
134
135
136	write(stringwftime,'(I7.7)')abs(wftime(indj))
137
138	idiagaa=0 !diagnostic on reading and opening files if 1 write more info
139	if (idiagaa.eq.1) then
140	!open(71,file='..\windtrans\data_type.txt') ! for testing: mc
141	!open(72,file='..\windtrans\list_windfield.txt') ! for testing: mc
142	open(unitdiagnostic1,file='list_windfield.txt') ! for testing: mc
143	!open(73,file='..\windtrans\list_global_att.txt') ! for testing: mc
144	open(unitdiagnostic2,file='list_global_att.txt') ! for testing: mc
145	!open(74,file='..\windtrans\list_variable_att.txt') ! for testing: mc
146	open(unitdiagnostic3,file='list_variable_att.txt') ! for testing: mc
147	!open(75,file='..\windtrans\seq_diagnostict.txt') ! for testing: mc
148	!open(75,file='..\windtrans\seq_diagnostict.txt') ! for testing: mc
149	!open(76,file='..\windtrans\vertical_wind_nov2_prima'//stringwftime//'.dat') ! for testing: mc
150	!open(77,file='..\windtrans\vertical_wind_nov2_dopo'//stringwftime//'.dat') ! for testing: mc
151	end if
152
153	gotgrid=0
154	ierr=0
155
156	!c-------------- open 4D (3D plus dtime dimension) wind meteo file
157	fnamenc=path(3)(1:length(3))&
158	//trim(wfname(indj))
159
160	iret = nf_open( fnamenc, NF_NOWRITE, ncid )
161	if (iret .ne. nf_noerr) then
162	write(*,9100) 'error doing open', fnamenc
163	ierr = -1
164	stop
165	end if
166	!c-------------- get information on dimension
167	iret = nf_inq( ncid, &
168	ndims_tot, nvars_tot, natts_tot, idimid_unlim )
169	if (iret .ne. nf_noerr) then
170	write(*,9100) 'error inquiring dimensions', fnamenc
171	ierr = -2
172	stop
173	end if
174
175	!c-------------inquiring abouot dimensions name -------
176	dimtimenum=0
177	do i = 1, min(ndims_tot,maxdim)
178	iret = nf_inq_dim( ncid, i, dimname(i), lendim(i) )
179	if (iret .ne. nf_noerr) then
180	write(*,9110) 'error inquiring dimensions for dim#', i, fnamenc
181	ierr = -2
182	stop
183	end if
184	if (dimname(i).eq.'time') dimtimenum=i
185	end do
186	!c-----------------------------------------------------------
187	!C------------inquiring about global attributes ---------
188
189	if (idiagaa .gt. 0) then
190	write(unitdiagnostic2,*)
191	write(unitdiagnostic2,*) 'attribute #, name, type, value'
192	end if
193	do 3401 iatt = 1, natts_tot
194	iret = nf_inq_attname( ncid, nf_global, iatt, attname)
195	if (iret .ne. nf_noerr) goto 3601
196	iret = nf_inq_att( ncid, nf_global, attname, ivtype, lenatt )
197	if (iret .ne. nf_noerr) goto 3601
198	if (ivtype .eq. 2) then
199	iret = nf_get_att_text( ncid, nf_global, attname, dumch1000 )
200	if (iret .ne. nf_noerr) goto 3601
201	i = max(1,min(1000,lenatt))
202
203	if (idiagaa .gt. 0) write(unitdiagnostic2,91010) &
204	iatt, attname(1:40), ivtype, lenatt, dumch1000(1:i)
205	else if (ivtype .eq. 4) then
206	iret = nf_get_att_int( ncid, nf_global, attname, idum )
207	if (iret .ne. nf_noerr) goto 3601
208	if (idiagaa .gt. 0) write(unitdiagnostic2,91020) &
209	iatt, attname(1:40), ivtype, lenatt, idum
210	else if ((ivtype .eq. 5) .and. (lenatt .eq. 1)) then
211	iret = nf_get_att_real( ncid, nf_global, attname, duma )
212	if (iret .ne. nf_noerr) goto 3601
213	if (idiagaa .gt. 0) write(unitdiagnostic2,91030) &
214	iatt, attname(1:40), ivtype, lenatt, duma
215	else if ((ivtype .eq. 5) .and. (lenatt .gt. 1)) then
216	allocate( duma_alloc(lenatt) )
217	iret = nf_get_att_real( ncid, nf_global, attname, duma_alloc )
218	if (iret .ne. nf_noerr) goto 3601
219	if (idiagaa .gt. 0) then
220	write(unitdiagnostic2,91010) iatt, attname(1:40), ivtype, lenatt
221	write(unitdiagnostic2,91040) (duma_alloc(i), i=1,lenatt)
222	end if
223
224	deallocate( duma_alloc )
225	else
226	if (idiagaa .gt. 0) write(unitdiagnostic2,'(i4,1x,a,2(1x,i6))') &
227	iatt, attname(1:40), ivtype, lenatt
228	goto 3401
229	endif
230
231	3401 continue
232	91010 format( i4, 1x, a, 2(1x,i6), 1x, a )
233	91020 format( i4, 1x, a, 2(1x,i6), 1x, i10 )
234	91030 format( i4, 1x, a, 2(1x,i6), 1x, 1pe12.4 )
235	91040 format(( 12x, 5(1pe12.4) ))
236
237	91050 format(a,1x,2(1x,i5), 12(1x,e15.8))
238	goto 3901
239
240	3601 write(*,9110) 'error inquiring global attribute', iatt, fnamenc
241	!write(73,9110) 'error inquiring global attribute', iatt, fnamenc
242	stop
243
244	3901 continue
245
246
247	!c----------------some checks on the varibles--------------------
248	do ii=1, nvars_tot
249	varid=ii
250	iret = nf_inq_var ( ncid, varid, varname, xtype, &
251	ndims, dimids, natts)
252	if (iret .ne. nf_noerr) then
253	write (,) 'error in nf_inq_var'
254	stop
255	end if
256
257	do iatt = 1, natts
258	iret = nf_inq_attname( ncid, varid, iatt, attname)
259	if (iret .ne. nf_noerr) goto 3602
260	iret = nf_inq_att( ncid, varid, attname, ivtype, lenatt )
261	if (iret .ne. nf_noerr) goto 3602
262	if (ivtype .eq. 2) then
263	dumch1000=''
264	iret = nf_get_att_text( ncid, varid, attname, dumch1000 )
265	if (iret .ne. nf_noerr) goto 3602
266	i = max(1,min(1000,lenatt))
267	if (attname.eq.'units') then
268	units(ii)=dumch1000
269	end if
270
271	if (idiagaa .gt. 0) write(unitdiagnostic3,91010) &
272	iatt, attname(1:40), ivtype, lenatt, dumch1000(1:i)
273	else if (ivtype .eq. 4) then
274	iret = nf_get_att_int( ncid, varid, attname, idum )
275	if (iret .ne. nf_noerr) goto 3602
276	if (idiagaa .gt. 0) write(unitdiagnostic3,91020) &
277	iatt, attname(1:40), ivtype, lenatt, idum
278	else if ((ivtype .eq. 5) .and. (lenatt .eq. 1)) then
279	iret = nf_get_att_real( ncid, varid, attname, duma )
280	if (iret .ne. nf_noerr) goto 3602
281	if (idiagaa .gt. 0) write(unitdiagnostic3,91030) &
282	iatt, attname(1:40), ivtype, lenatt, duma
283	else if ((ivtype .eq. 5) .and. (lenatt .gt. 1)) then
284	allocate( duma_alloc(lenatt) )
285	iret = nf_get_att_real( ncid, varid, attname, duma_alloc )
286	if (iret .ne. nf_noerr) goto 3602
287	if (idiagaa .gt. 0) then
288	write(unitdiagnostic3,91010) iatt, attname(1:40), ivtype, lenatt
289	write(unitdiagnostic3,91040) (duma_alloc(i), i=1,lenatt)
290	end if
291
292	deallocate( duma_alloc )
293	else
294	if (idiagaa .gt. 0) write(unitdiagnostic3,'(i4,1x,a,2(1x,i6))') &
295	iatt, attname(1:40), ivtype, lenatt
296	exit
297	endif
298
299	end do
300	goto 3902
301	3602 write(*,9110) 'error inquiring variable attribute', varname, &
302	iatt, attname, fnamenc
303	!write(74,9110) 'error inquiring variable attribute', varname, &
304	!iatt, attname, fnamenc
305	stop
306	3902 continue
307
308	!write(75,*)'start',ncid, varid, varname, xtype, & !more diagnostic for testing mc
309	!ndims, dimids, natts,'fine'
310
311	xtypev(ii)=xtype
312	varnamev(ii)=varname
313	ndimsv(ii)=ndims
314	do j=1,maxdim
315	dimidsm(ii,j)=dimids(j)
316	end do
317	!write(71,9031)i,varnamev(ii),xtypev(ii),ndims,units(ii) !WRITE diagnostic to file BY M C
318
319
320
321	end do
322
323	!C-------------- find the time --------------
324
325	varname='date' !read date
326	call check_variable(varname,fnamenc,maxdim,nf_int, &
327	id_var,ndims,id_dim,ierr,ncid)
328	if (ierr.lt.0) goto 100
329	lendim_exp=0
330	if (ndims.eq.0) then !this shoudl happen if there is only one time per file
331	ndims=1
332	lendim_exp(1)=1
333	else
334	do j = 1, ndims
335	lendim_exp(j) = lendim(id_dim(j))
336	end do
337	end if
338	call allocatedumarray(ndims,lendim_exp,maxdim,nf_int)
339	do j=1, ndims
340	istart(j) = 1 !ndims
341	icount(j) = lendim_exp(j)
342	end do
343	!here we assume taht ndims=1 for date variable
344	iret = &
345	nf_get_vara_int( ncid, id_var, istart, icount, dumarray1D_int)
346	if (iret .ne. nf_noerr) then
347	write(*,9100) 'error inquiring var value ' // varname, fnamenc
348	ierr = -5
349	goto 100
350	end if
351	do i=1,lendim_exp(1) !one dimension only expected!
352	date_aid(i)=dumarray1D_int(i)
353	end do
354
355
356	varname='datesec' !read date in seconds
357	call check_variable(varname,fnamenc,maxdim,nf_int, &
358	id_var,ndims,id_dim,ierr,ncid)
359	if (ierr.lt.0) goto 100
360	lendim_exp=0
361	if (ndims.eq.0) then
362	ndims=1
363	lendim_exp(1)=1
364	else
365	do j = 1, ndims
366	lendim_exp(j) = lendim(id_dim(j))
367	end do
368	end if
369	call allocatedumarray(ndims,lendim_exp,maxdim,nf_int)
370	do j=1, ndims
371	istart(j) = 1 !ndims
372	icount(j) = lendim_exp(j)
373	end do
374	iret = &
375	nf_get_vara_int( ncid, id_var, istart, icount, dumarray1D_int)
376	if (iret .ne. nf_noerr) then
377	write(*,9100) 'error inquiring var value ' // varname, fnamenc
378	ierr = -5
379	goto 100
380	end if
381
382	do i=1,lendim_exp(1) !one dimension only expected
383	jul=juldate(date_aid(i),0)
384
385	jul=jul+real(dumarray1D_int(i),kind=dp)/86400._dp
386	dumb= (jul-bdate)*86400._dp
387	idumb=nint(dumb)
388	if (idumb.eq.wftime(indj)) exit
389	end do
390	itime=i !until here works.... it select right file and record....
391	!time_interval=wftime(indj)-wftime(indj-1)
392	!c-------- finish find time --------------
393	if(idiagaa.eq.1) then
394	write (unitdiagnostic1,*)wfname(indj),wftime(indj),itime,n !write diagnostic to file for test reason by, by mc
395	end if
396
397	!c---------------- Inquiring about varnames and real fields ---------------------
398	strswitch=.false. ! shoud become .true. below if wind file has the right data
399	hflswitch=.false. ! shoud become .true. below if wind file has the right data
400	nvar_exp_in_meteo_field=17
401	varnamev(1)='U' !horizontal wind m/s
402	varnamev(2)='V' !horizontal wind m/s
403	varnamev(3)='OMEGA' ! vertical wind Pa/s
404	varnamev(4)='T' !temperature k
405	varnamev(5)='Q' ! Specific humidity kg/kg
406	varnamev(6)='PS' !surface pressure Pa
407	varnamev(7)='CLDTOT' !total cloud cover fraction
408	varnamev(8)='TREFHT' !2m temperature k
409	varnamev(9)='PRECL' !large scale precipitation m/s must become mm/hopur in FLEXPARTECMWF ecmwf 142 but deaccumulated mm/hour /I must multiply by output time in s dvivide by hours and multipluy by 1000
410	varnamev(10)='PRECC' !convective precipitation ms/ must become mm/hour in FLEXPARTECMWF cmwf 143 but deaccumulated mm/hour /I must multiply by output time in s dvivide by hours and multipluy by 1000
411	varnamev(11)='SHFLX' !sensible heat fluxes W/m^2 in FLEXPARTECMWF ecmwf 146 but deaccumulated w/m^2 / this is already in Watt m^2 so nothing to do
412	varnamev(12)='TAUX' !surface stress east-west N/m^2 in FLEXPARTECMWF ecmwf 180 but deaccumulated to N/m^2 / this is already in N/m^2 so nothign to do
413	varnamev(13)='TAUY' !surface stress north-south N/m^2 in FLEXPARTECMWF ecmwf 181 but deaccumulated to N/m^2 this is already in N/m^2 so nothing to do
414	varnamev(14)='U10' !ten meters wind speed m/s
415	varnamev(15)='QREFHT' !speciifc humidity
416	varnamev(16)='SNOWHLND' !m water equivakent snow depth !double check if must be also summed SNOWHICE
417	varnamev(17)='FSDS' !downwelling solar flux atsurface! double checck if this correct for stomata opeining parameterizations !note that in ecmwf SSR ecmwf 176 but deaccumulated
418
419	!---------------- check and load file contents
420	do i=1,nvar_exp_in_meteo_field !variable_loop
421	varname=varnamev(i)
422	call check_variable(varname,fnamenc,maxdim,nf_float, &
423	id_var,ndims,id_dim,ierr,ncid)
424	if (ierr.ne.0) goto 100
425	lendim_exp=0
426	indextime=0
427	do j = 1, ndims
428	lendim_exp(j) = lendim(id_dim(j))
429	if (id_dim(j).eq.dimtimenum) indextime=j !this indextime select whcih is the dimension associated with time
430	end do
431	if (dimtimenum.eq.0.and.indextime.eq.0) then !this means that there is no time dimension in the file and therefore we create an extra time dimension of leght 1.
432	!notimedim=1
433	ndims=ndims+1
434	lendim_exp(ndims)=1
435	indextime=ndims
436	end if
437
438	if (indextime.ne.ndims) stop 'ERROR readwind the time is expected &
439	to be always the last dimension in the stored variables'
440
441	call allocatedumarray(ndims,lendim_exp,maxdim,nf_float)
442
443	do j=1, ndims
444	istart(j) = 1 !ndims
445	icount(j) = lendim_exp(j)
446	end do
447
448	if (ndims.eq.1) then
449	iret = &
450	nf_get_vara_real( ncid, id_var, istart, icount, dumarray1D_real)
451	if (iret .ne. nf_noerr) then
452	write(*,9100) 'error inquiring var value ' // varname, fnamenc
453	ierr = -5
454	goto 100
455	end if
456	else if (ndims.eq.2) then
457	iret = &
458	nf_get_vara_real( ncid, id_var, istart, icount, dumarray2D_real)
459	if (iret .ne. nf_noerr) then
460	write(*,9100) 'error inquiring var value ' // varname, fnamenc
461	ierr = -5
462	goto 100
463	end if
464	else if (ndims.eq.3) then
465	iret = &
466	nf_get_vara_real( ncid, id_var, istart, icount, dumarray3D_real)
467	if (iret .ne. nf_noerr) then
468	write(*,9100) 'error inquiring var value ' // varname, fnamenc
469	ierr = -5
470	if ((varname.eq.'TAUx').or.(varname.eq.'TAUY').or.(varname.eq.&
471	'SHFLX')) then
472	print ,'***********************************'
473	print ,' no stress available **************'
474	print ,'***********************************'
475	cycle ! skyp the rest and go to the next variable
476	end if
477	goto 100
478	end if
479	else if (ndims.eq.4) then
480	iret = &
481	nf_get_vara_real( ncid, id_var, istart, icount, dumarray4D_real)
482	if (iret .ne. nf_noerr) then
483	write(*,9100) 'error inquiring var value ' // varname, fnamenc
484	ierr = -5
485	goto 100
486	end if
487	end if
488	!write(75,*) varname
489
490	do jy=0,nymin1
491	do ix=0, nxfield-1
492
493	if (varname.eq.'U') then
494	do kz=1,nuvz-1 !recall that there is an augmented surface layer where U10 will be stored
495	uuh(ix,jy,nuvz+1-kz)=dumarray4D_real(ix+1,jy+1,kz,itime)
496	end do
497	else if (varname.eq.'V') then
498	do kz=1,nuvz-1 !recall taht there is an augmented surface layer where V10 will be stored
499	vvh(ix,jy,nuvz+1-kz)=dumarray4D_real(ix+1,jy+1,kz,itime)
500	end do
501	else if (varname.eq.'OMEGA') then
502	do kz=1,nuvz-1 !recall that there is an augmented surface layer where W must be defined
503	wwh(ix,jy,nuvz+1-kz)=dumarray4D_real(ix+1,jy+1,kz,itime)
504	end do
505	else if (varname.eq.'T') then
506	do kz=1,nuvz-1 !recall that there is an augmented surface layer where T must be defined
507	tth(ix,jy,nuvz+1-kz,n)=dumarray4D_real(ix+1,jy+1,kz,itime)
508	end do
509	else if (varname.eq.'Q') then
510	do kz=1,nuvz-1 !recall that there is an augmented surface layer where Q must be defined
511	qvh(ix,jy,nuvz+1-kz,n)=dumarray4D_real(ix+1,jy+1,kz,itime)
512	if (qvh(ix,jy,nuvz+1-kz,n).lt.0) qvh(ix,jy,nuvz+1-kz,n)=0.
513	end do
514	else if (varname.eq.'PS') then
515	ps(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)
516	else if (varname.eq.'CLDTOT') then
517	tcc(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)
518	else if (varname.eq.'TREFHT') then
519	tt2(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)
520	else if (varname.eq.'PRECL') then
521	lsprec(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)* &
522	3600.*1000. !convet from m/s to mm/hour this must be controlled and eventually replaced with cumulated values
523	else if (varname.eq.'PRECC') then
524	convprec(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)* &
525	3600.*1000. !convet from m/s to mm/hour
526	else if (varname.eq.'SHFLX') then
527	sshf(ix,jy,1,n)=-dumarray3D_real(ix+1,jy+1,itime) ! note the minus sign here because FLEXPART assume positive FLUX downward like ECMWF: see e.g.http://tigge.ecmwf.int/tigge/d/show_object/table=parameters/name=time_integrated_surface_sensible_heat_flux/levtype=sfc/
528	hflswitch=.true. ! Heat flux available
529	else if (varname.eq.'TAUX') then
530	ewss(ix,jy)=dumarray3D_real(ix+1,jy+1,itime)
531	strswitch=.true. ! put this to true if stress are avilable
532	else if (varname.eq.'TAUY') then
533	nsss(ix,jy)=dumarray3D_real(ix+1,jy+1,itime)
534	else if (varname.eq.'U10') then ! in NORESM U10 contain sped not component so it is partitioned using ground stresses, see below
535	if (strswitch) then
536	u10(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)* ewss(ix,jy)/ &
537	sqrt(ewss(ix,jy)2+nsss(ix,jy)2)
538	v10(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)* nsss(ix,jy)/ &
539	sqrt(ewss(ix,jy)2+nsss(ix,jy)2)
540	else
541	u10(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime) ! if stress not available then see below
542	v10(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)
543	end if
544	else if (varname.eq.'QREFHT') then !
545	td2(ix,jy,1,n)=dewpoint(tt2(ix,jy,1,n), &
546	dumarray3D_real(ix+1,jy+1,itime),ps(ix,jy,1,n))
547	qv2(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime)
548	!else if (varname.eq.'SOLARRADIATION') then !
549	! ssr(ix,jy,1,n)=
550	else if (varname.eq.'SNOWHLND') then
551	sd(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime) !m of water equivalent
552	else if (varname.eq.'FSDS') then
553	ssr(ix,jy,1,n)=dumarray3D_real(ix+1,jy+1,itime) !downwelling solar radition flux in w/m^2
554	end if
555
556	end do
557	end do
558
559	!c----------- deallocation added my mc on 12-2013 --------
560	if (allocated(dumarray4D)) deallocate(dumarray4D)
561	if (allocated(dumarray3D)) deallocate(dumarray3D)
562	if (allocated(dumarray2D)) deallocate(dumarray2D)
563	if (allocated(dumarray1D)) deallocate(dumarray1D)
564	if (allocated(dumarray4D_real)) deallocate(dumarray4D_real)
565	if (allocated(dumarray3D_real)) deallocate(dumarray3D_real)
566	if (allocated(dumarray2D_real)) deallocate(dumarray2D_real)
567	if (allocated(dumarray1D_real)) deallocate(dumarray1D_real)
568	if (allocated(dumarray4D_int)) deallocate(dumarray4D_int)
569	if (allocated(dumarray3D_int)) deallocate(dumarray3D_int)
570	if (allocated(dumarray2D_int)) deallocate(dumarray2D_int)
571	if (allocated(dumarray1D_int)) deallocate(dumarray1D_int)
572
573	!c-------------------------------------------------------------------------
574	!c------------- read number of vertical level
575	!c------------- note taht U, V, T, Q, OMEGA are colocated in CAM3.0/CAM4.0 see user's guide to thE NCAR CAM 3.0 page 38-
576	!c------------- so nwz and nuvz are teh same here while tadot levels form surface to top are nwz+1=27 for CAm 4.+0!
577	if (varname.eq.'OMEGA') then
578	if (nwz.ne.lendim_exp(3)) stop 'READWIND nwz NOT CONSISTENT'
579	else if (varname.eq.'U') then
580	if (nuvz-1.ne.lendim_exp(3)) stop 'READWIND nwz NOT CONSITENT'
581	end if
582	!c----------
583
584	end do !variable_loop !
585
586	!c---------------------------
587	iret = nf_close( ncid ) !closee netcdf files
588
589	!c----------------------- more diagnostic for vertical velcoity ----------
590	! do jy=0,nymin1
591	! do ix=0, nxfield-1
592	! do kz=1,26
593	! write (76,'(3f12.4,4E15.6)')0+dxix,ylat0+dyjy,kz*1.,wwh(ix,jy,kz)
594	! end do
595	! end do
596	! end do
597
598	!c----------------- here we have to put the call to the routine that will transform OMEGA in etadot * delta_p/delta_eta
599	!deltat=wftime(indj+1)-wftime(indj-1) !time interval around indj
600	if (method_w.eq.1) then ! <<<<---------------------if on method_w
601	call transform_omega_etadot(uuh,vvh,wwh,n)
602	end if
603	!C-----C more dignostic for vertical velocity ---------------
604	!do jy=0,nymin1
605	! do ix=0, nxfield-1
606	! do kz=1,26
607	! write (77,'(3f12.4,4E15.6)')0+dxix,ylat0+dyjy,kz*1.,wwh(ix,jy,kz)
608	! end do
609	! end do
610	!end do
611	!c-------------------------------------------------------------
612
613	!c------------------------------------------------------------------------------------------------------------------
614
615
616
617	!C---- If desired, shift all grids by nxshift grid cells ------
618
619	if (xglobal) then
620	call shift_field_0(ewss,nxfield,ny)
621	call shift_field_0(nsss,nxfield,ny)
622	!call shift_field_0(oro,nxfield,ny)
623	!call shift_field_0(excessoro,nxfield,ny)
624	!call shift_field_0(lsm,nxfield,ny)
625	call shift_field(ps,nxfield,ny,1,1,2,n)
626	call shift_field(sd,nxfield,ny,1,1,2,n)
627	call shift_field(msl,nxfield,ny,1,1,2,n)
628	call shift_field(tcc,nxfield,ny,1,1,2,n)
629	call shift_field(u10,nxfield,ny,1,1,2,n)
630	call shift_field(v10,nxfield,ny,1,1,2,n)
631	call shift_field(tt2,nxfield,ny,1,1,2,n)
632	call shift_field(td2,nxfield,ny,1,1,2,n)
633	call shift_field(lsprec,nxfield,ny,1,1,2,n)
634	call shift_field(convprec,nxfield,ny,1,1,2,n)
635	call shift_field(sshf,nxfield,ny,1,1,2,n)
636	call shift_field(ssr,nxfield,ny,1,1,2,n)
637	call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
638	call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
639	call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
640	call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
641	call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
642	endif
643	!C**********************************************
644
645	do i=0,nxmin1
646	do j=0,nymin1
647	surfstr(i,j,1,n)=sqrt(ewss(i,j)2+nsss(i,j)2)
648	end do
649	end do
650
651	!c*****************************************************
652	if ((.not.hflswitch).or.(.not.strswitch)) then
653	write(,) 'WARNING: No flux data contained in wind file ', &
654	wfname(indj)
655
656	! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
657	! As ECMWF has increased the model resolution, such that now the first model
658	! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
659	! (3rd model level in FLEXPART) for the profile method
660	!***************************************************************************
661
662	do i=0,nxmin1
663	do j=0,nymin1
664	plev1=akz(3)+bkz(3)*ps(i,j,1,n)
665	pmean=0.5*(ps(i,j,1,n)+plev1)
666	tv=tth(i,j,3,n)(1.+0.61qvh(i,j,3,n))
667	fu=-r_air*tv/ga/pmean
668	hlev1=fu*(plev1-ps(i,j,1,n)) ! HEIGTH OF FIRST MODEL LAYER
669	ff10m= u10(i,j,1,n) !note that in NORESM, u10 and v10, if strswitch.eq.false, contain the speed not the components
670	fflev1=sqrt(uuh(i,j,3)2+vvh(i,j,3)2)
671	call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
672	tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
673	surfstr(i,j,1,n),sshf(i,j,1,n))
674	if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
675	if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
676	end do
677	end do
678	do i=0,nxmin1
679	do j=0,nymin1
680	u10(i,j,1,n)=u10(i,j,1,n)uuh(i,j,2)/ & ! *NOTE:since we had no taux and tauy we used speed at 10 meters to determine total stress see above
681	sqrt(uuh(i,j,2)2+vvh(i,j,2)2) ! **therefore to assgin the wind component at 10 meters we have to do a further assumption
682	v10(i,j,1,n)=v10(i,j,1,n)vvh(i,j,2)/ & ! *this assumtpion is to take direction equal to the one of the first computed wind level
683	sqrt(uuh(i,j,2)2+vvh(i,j,2)2) ! **which in general is probably not that correct..
684	end do
685	end do
686
687	endif ! on surface stresses
688
689	! Assign 10 m wind to model level at eta=1.0 to have one additional model
690	! level at the ground
691	! Specific humidity is taken the same as at one level above
692	! Temperature is taken as 2 m temperature
693	!**************************************************************************
694
695	do i=0,nxmin1
696	do j=0,nymin1
697	wwh(i,j,1)=0. !NORESM in ECMWf the OMEGA (hybrid system velocity in Pa/s) value at the ground is used alternatively we shoudl use etadot (1/s) at the ground to obtain wwh at the ground
698	uuh(i,j,1)=u10(i,j,1,n)
699	vvh(i,j,1)=v10(i,j,1,n)
700	!qvh(i,j,1,n)=qvh(i,j,2,n) !ECMWF version
701	qvh(i,j,1,n)=qv2(i,j,1,n) ! this replace the above in NORESM version
702	tth(i,j,1,n)=tt2(i,j,1,n)
703	end do
704	end do
705
706
707	nlev_ec=nuvz-1 !in ECMWF version nelev_ec is the counter of ECMWF number of vertical levels for U (e.g. 91 with 92 levels)
708
709	return
710
711
712
713	100 continue
714	print *,'error reading',varname,'erros code',ierr
715
716	stop
717
718
719	return
720	end
721
722

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

Download in other formats:

Original Format