Context Navigation

gridcheck.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: 29.8 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	Subroutine gridcheck()
26
27	use noresm_variables
28	use par_mod
29	use com_mod
30	use conv_mod
31	use cmapf_mod, only: stlmbr,stcm2p
32
33	implicit none
34
35	include 'netcdf.inc'
36
37	!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
38	!* *
39	!* Read grid strucutre from netcdf files *
40	!* set grid for FLEXPART-NorESM *
41	!* check if variables are there in the netcdf files *
42	!* note: part of netcdf reading strucutre adapted from the routines by *
43	!* Fast and Easter in FLEXPART-WRF *
44	!* *
45	!* *
46	!* Author: *
47	!* M. Cassiani 2016 *
48	!* *
49	!* *
50	!* *
51	!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
52
53
54	!integer maxdim,maxvar,maxtime !moved in par_mod.f90
55	real(kind=4) eps
56
57	!parameter(maxdim=9,maxvar=70,maxtime=12, !this must be moved in par_mod when finalized
58	parameter(eps=0.0001)
59
60	integer :: gotGrid
61	real(kind=4) :: sizesouth,sizenorth,xauxa,pint
62
63	!c---------------------------------------------
64	real(kind=4) :: xaux1,xaux2,yaux1,yaux2
65	real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
66	integer :: nyfieldin,nxfieldin
67
68	!C-------------- declaration for netcdf reading
69
70	real(kind=4) :: duma
71	real(kind=4), allocatable, dimension(:) :: duma_alloc
72
73	integer :: ierr !error code message
74	integer :: idiagaa !flag
75	integer :: id_var,id_dim(maxdim)
76	integer :: nvar_exp_in_grid_atm_nc,nvar_exp_in_meteo_field
77
78	integer :: i,j,iatt,idimid_unlim,idum,iret,ivtype
79	integer :: ix,jy,ifn
80	integer :: lenatt,lendim(maxdim)
81	integer :: natts_tot,ncid,ndims_tot,nvars_tot
82	integer :: sizetype
83	character*110 :: fnamenc
84	character*80 :: dimname(maxdim)
85	character*80 :: attname
86	character*160 :: varname,vartype
87	character*160 :: varnamev(maxvar)
88	character*160 :: vartypev(maxvar)
89
90	integer :: ndimsv(maxvar)
91	integer :: dimidsv(maxvar,maxdim)
92	character*1000 :: dumch1000
93
94	!integer istart(maxdim),icount(maxdim)
95	integer :: xtype,xtypev(maxvar)
96	integer :: ndims
97	integer :: dimids(maxdim),dimidsm(maxvar,maxdim)
98	integer :: LENDIM_EXP(maxdim),LENDIM_MAX(maxdim)
99	integer :: natts
100	integer :: varid
101
102	!*************************************************************************************************
103	if(ideltas.gt.0) then
104	ifn=1
105	else
106	ifn=numbwf
107	endif
108
109	!*************************************************************************************************
110	9100 format( / '** read_noresmout_gridinfo -- ', a / &
111	'file = ', a )
112	9110 format( / '** read_noresmout_gridinfo -- ', a, 1x, i8 / &
113	'file = ', a )
114	9120 format( / '** read_noresmout_gridinfo -- ', a, 2(1x,i8) / &
115	'file = ', a )
116
117	9030 format( a, 2i6, 2(2x,a) )
118	9031 format( 1i4,1a20,2i4)
119
120
121	idiagaa=0 !diagnostic on reading and opening files if 1 write more info
122	if (idiagaa.eq.1) then
123	!!open(71,file='..\options\data_type.txt') ! for testing: mc
124	!open(73,file='..\options\list_global_att.txt') ! for testing: mc
125	open(unitdiagnostic2,file='list_global_att.txt') ! for testing: mc
126	!!open(75,file='..\options\seq_diagnostict.txt') ! for testing: mc
127	end if
128
129	gotgrid=0
130	ierr=0
131	!c-------------- open grid structure file
132
133	fnamenc=path(5)(1:length(5)) !'..\..\..\data_set\grid_atm.nc' !grid_atm.nc'
134	ncid = 10
135	iret = nf_open( fnamenc, NF_NOWRITE, ncid )
136	if (iret .ne. nf_noerr) then
137	write(*,9100) 'error doing open', fnamenc
138	ierr = -1
139	end if
140	!c-------------- get information on dimension
141
142	iret = nf_inq( ncid, &
143	ndims_tot, nvars_tot, natts_tot, idimid_unlim )
144	if (iret .ne. nf_noerr) then
145	write(*,9100) 'error inquiring dimensions', fnamenc
146	ierr = -2
147	stop
148	end if
149
150
151	!c-------------inquiring abouot dimensions name -------
152	do i = 1, min(ndims_tot,maxdim)
153	iret = nf_inq_dim( ncid, i, dimname(i), lendim(i) )
154	if (iret .ne. nf_noerr) then
155	write(*,9110) 'error inquiring dimensions for dim#', i, fnamenc
156	ierr = -2
157	stop
158	end if
159	end do
160
161	!c------------inquiring about global attributes ---------
162
163	if (idiagaa .gt. 0) then
164	write(73,*)
165	write(73,*) 'grid strucuture: attribute #, name, type, value'
166	end if
167	do 3400 iatt = 1, natts_tot
168	iret = nf_inq_attname( ncid, nf_global, iatt, attname)
169	if (iret .ne. nf_noerr) goto 3600
170	iret = nf_inq_att( ncid, nf_global, attname, ivtype, lenatt )
171	if (iret .ne. nf_noerr) goto 3600
172	if (ivtype .eq. 2) then
173	iret = nf_get_att_text( ncid, nf_global, attname, dumch1000 )
174	if (iret .ne. nf_noerr) goto 3600
175	i = max(1,min(1000,lenatt))
176	if (idiagaa .gt. 0) write(73,91010) &
177	iatt, attname(1:40), ivtype, lenatt, dumch1000(1:i)
178	else if (ivtype .eq. 4) then
179	iret = nf_get_att_int( ncid, nf_global, attname, idum )
180	if (iret .ne. nf_noerr) goto 3600
181	if (idiagaa .gt. 0) write(73,91020) &
182	iatt, attname(1:40), ivtype, lenatt, idum
183	else if ((ivtype .eq. 5) .and. (lenatt .eq. 1)) then
184	iret = nf_get_att_real( ncid, nf_global, attname, duma )
185	if (iret .ne. nf_noerr) goto 3600
186	if (idiagaa .gt. 0) write(73,91030) &
187	iatt, attname(1:40), ivtype, lenatt, duma
188	else if ((ivtype .eq. 5) .and. (lenatt .gt. 1)) then
189	allocate( duma_alloc(lenatt) )
190	iret = nf_get_att_real( ncid, nf_global, attname, duma_alloc )
191	if (iret .ne. nf_noerr) goto 3600
192	if (idiagaa .gt. 0) then
193	write(73,91010) iatt, attname(1:40), ivtype, lenatt
194	write(73,91040) (duma_alloc(i), i=1,lenatt)
195	end if
196	deallocate( duma_alloc )
197	else
198	if (idiagaa .gt. 0) write(73,'(i4,1x,a,2(1x,i6))') &
199	iatt, attname(1:40), ivtype, lenatt
200	goto 3400
201	endif
202
203	3400 continue
204	91010 format( i4, 1x, a, 2(1x,i6), 1x, a )
205	91020 format( i4, 1x, a, 2(1x,i6), 1x, i10 )
206	91030 format( i4, 1x, a, 2(1x,i6), 1x, 1pe12.4 )
207	91040 format(( 12x, 5(1pe12.4) ))
208
209	91050 format(a,1x,2(1x,i5), 12(1x,e15.8))
210
211	goto 3900
212
213	3600 write(*,9110) 'error inquiring attribute', iatt, fnamenc
214	write(73,9110) 'error inquiring attribute', iatt, fnamenc
215	stop
216
217	3900 continue
218
219
220	!c---------------- Inquiring about varnames and real fields ---------------------
221	do i=1, nvars_tot
222	varid=i
223	iret = nf_inq_var ( ncid, varid, varname, xtype, &
224	ndims, dimids, natts)
225	if (iret .ne. nf_noerr) then
226	write (,) 'error in nf_inq_var'
227	stop
228	end if
229	!if (idiagaa.eq.1) then
230	!!write(75,*)'grid structure inquire',ncid, varid, varname, xtype, &
231	!!ndims, dimids, natts,'fine'
232	!end if
233
234	xtypev(i)=xtype
235	varnamev(i)=varname
236	ndimsv(i)=ndims
237	do j=1,maxdim
238	dimidsm(i,j)=dimids(j)
239	end do
240	!if (idiagaa.eq.1) then
241	!!write(71,9031)i,varnamev(i),xtypev(i),ndims !WRITE diagnostic to file
242	!end if
243	end do
244	! stop
245
246	!c---------- number of varibles to recover from netcdf file
247	nvar_exp_in_grid_atm_nc=6
248	!c--------- 1D fields
249	varnamev(1)='lat'
250	varnamev(2)='lon'
251	!C--------- 2D fields
252	varnamev(3)='PHIS'
253	varnamev(4)='LANDFRAC'
254	varnamev(5)='SGH'
255	varnamev(6)='xv'
256
257	!c--------- read fields & set grid --------------------------------
258	do i=1,nvar_exp_in_grid_atm_nc
259	varname=varnamev(i)
260	call check_variable(varname,fnamenc,maxdim,nf_double, &
261	id_var,ndims,id_dim,ierr,ncid)
262	if (ierr.lt.0) goto 100
263	lendim_exp=0
264	do j = 1, ndims
265	lendim_exp(j) = lendim(id_dim(j))
266	end do
267
268	call allocatedumarray(ndims,lendim_exp,maxdim,nf_double)
269	do j=1, ndims
270	istart(j) = 1 !ndims
271	icount(j) = lendim_exp(j)
272	end do
273	if (ndims.eq.1) then
274	iret = &
275	nf_get_vara_double( ncid, id_var, istart, icount, dumarray1D)
276	if (iret .ne. nf_noerr) then
277	write(*,9100) 'error inquiring var value ' // varname, fnamenc
278	ierr = -5
279	goto 100
280	end if
281	else if (ndims.eq.2) then
282	iret = &
283	nf_get_vara_double( ncid, id_var, istart, icount, dumarray2D)
284	if (iret .ne. nf_noerr) then
285	write(*,9100) 'error inquiring var value ' // varname, fnamenc
286	ierr = -5
287	goto 100
288	end if
289	else if (ndims.eq.3) then
290	iret = &
291	nf_get_vara_double( ncid, id_var, istart, icount, dumarray3D)
292	if (iret .ne. nf_noerr) then
293	write(*,9100) 'error inquiring var value ' // varname, fnamenc
294	ierr = -5
295	goto 100
296	end if
297	end if
298
299	!if (idiagaa.eq.1) then
300	!!write(75,*)'gridcheck',varname
301	!end if
302
303	!here load the arrays
304	if (varname.eq.'lat') then
305	nyfieldin=lendim_exp(1)
306	yaux2in=dumarray1D(lendim_exp(1)) !last point
307	yaux1in=dumarray1D(1) !first point
308	ny=nyfieldin
309	else if (varname.eq.'lon') then
310	nxfieldin=lendim_exp(1)
311	xaux2in=dumarray1D(lendim_exp(1)) !last point
312	xaux1in=dumarray1D(1) !first point
313	nxfield=nxfieldin
314
315	!c-------------- now set grid structure --------------
316	xaux1=xaux1in
317	xaux2=xaux2in
318	yaux1=yaux1in
319	yaux2=yaux2in
320
321
322	if (xaux1.gt.180.) xaux1=xaux1-360.0
323	if (xaux2.gt.180.) xaux2=xaux2-360.0
324	if (xaux1.lt.-180.) xaux1=xaux1+360.0
325	if (xaux2.lt.-180.) xaux2=xaux2+360.0
326	if (xaux2.lt.xaux1) xaux2=xaux2+360.0
327	xlon0=xaux1
328	ylat0=yaux1
329	!ny=nyfieldin
330	dx=(xaux2-xaux1)/real(nxfield-1)
331	dy=(yaux2-yaux1)/real(ny-1)
332	dxconst=180./(dxr_earthpi)
333	dyconst=180./(dyr_earthpi)
334	gotGrid=1
335
336	! Check whether fields are global
337	! If they contain the poles, specify polar stereographic map
338	! projections using the stlmbr- and stcm2p-calls
339	!***********************************************************
340
341	xauxa=abs(xaux2+dx-360.-xaux1)
342	if (xauxa.lt.0.001) then
343	nx=nxfield+1 ! field is cyclic
344	xglobal=.true.
345	if (abs(nxshift).ge.nx) &
346	stop 'nxshift in file par_mod is too large'
347	xlon0=xlon0+real(nxshift)*dx
348	else
349	nx=nxfield
350	xglobal=.false.
351	if (xglobal.eqv..false.) &
352	stop 'Noresm/CAM simualation supposed to have global wind field and &
353	must use nxshift different from zero'
354	endif
355	nxmin1=nx-1
356	nymin1=ny-1
357	if (xlon0.gt.180.) xlon0=xlon0-360.
358	xauxa=abs(yaux1+90.)
359	if (xglobal.and.xauxa.lt.0.001) then
360	sglobal=.true. ! field contains south pole
361	! Enhance the map scale by factor 3 (*2=6) compared to north-south
362	! map scale
363	sizesouth=6.*(switchsouth+90.)/dy
364	call stlmbr(southpolemap,-90.,0.)
365	call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
366	sizesouth,switchsouth,180.)
367	switchsouthg=(switchsouth-ylat0)/dy
368	else
369	sglobal=.false.
370	switchsouthg=999999.
371	endif
372	xauxa=abs(yaux2-90.)
373	if (xglobal.and.xauxa.lt.0.001) then
374	nglobal=.true. ! field contains north pole
375	! Enhance the map scale by factor 3 (*2=6) compared to north-south
376	! map scale
377	sizenorth=6.*(90.-switchnorth)/dy
378	call stlmbr(northpolemap,90.,0.)
379	call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
380	sizenorth,switchnorth,180.)
381	switchnorthg=(switchnorth-ylat0)/dy
382	else
383	nglobal=.false.
384	switchnorthg=999999.
385	endif
386	if (nxshift.lt.0) &
387	stop 'nxshift (par_mod) must not be negative'
388	if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
389
390	!c--------------- end of horizontal grid structure definitions --------
391
392
393	else if (varname.eq.'PHIS') then
394	do jy=0,ny-1
395	do ix=0,nxfield-1
396	oro(ix,jy)=sngl(dumarray2D(ix+1,jy+1))/ga
397	end do
398	end do
399	else if (varname.eq.'SGH') then
400	do jy=0,ny-1
401	do ix=0,nxfield-1
402	excessoro(ix,jy)=sngl(dumarray2D(ix+1,jy+1))
403	end do
404	end do
405	else if (varname.eq.'LANDFRAC') then
406	do jy=0,ny-1
407	do ix=0,nxfield-1
408	lsm(ix,jy)=sngl(dumarray2D(ix+1,jy+1))
409	end do
410	end do
411	end if
412
413	end do
414	!c------------------ set the grid ---------------------
415
416	print *,'gotgrid',gotgrid
417	!error message if no fields found with correct first longitude in it
418	if (gotGrid.eq.0) then
419	print,'**ERROR: horizontal grid strucutre not defined'
420	stop
421	endif
422
423	! goto 200
424	iret = nf_close( ncid )
425	!c-------------- open 4D (3D plus time dimension) wind meteo file
426
427
428	fnamenc=path(3)(1:length(3))//trim(wfname(ifn))
429
430	iret = nf_open( fnamenc, NF_NOWRITE, ncid )
431	if (iret .ne. nf_noerr) then
432	write(*,9100) 'error doing open', fnamenc
433	ierr = -1
434	!stop
435	end if
436	!c-------------- get information on dimension
437	iret = nf_inq( ncid, &
438	ndims_tot, nvars_tot, natts_tot, idimid_unlim )
439	if (iret .ne. nf_noerr) then
440	write(*,9100) 'error inquiring dimensions', fnamenc
441	ierr = -2
442	stop
443	end if
444	idiagaa=1
445	!c-------------inquiring abouot dimensions name -------
446	do i = 1, min(ndims_tot,maxdim)
447	iret = nf_inq_dim( ncid, i, dimname(i), lendim(i) )
448	if (iret .ne. nf_noerr) then
449	write(*,9110) 'error inquiring dimensions for dim#',i,fnamenc
450	ierr = -2
451	stop
452	end if
453	end do
454	!c-----------------------------------------------------------
455	!C------------inquiring about global attributes ---------
456
457	if (idiagaa .gt. 0) then
458	write(,)
459	write(73,*)'gridcheck-windfiled: attribute #, name, type, value'
460	end if
461	do 3401 iatt = 1, natts_tot
462	iret = nf_inq_attname( ncid, nf_global, iatt, attname)
463	if (iret .ne. nf_noerr) goto 3601
464	iret = nf_inq_att( ncid, nf_global, attname, ivtype, lenatt )
465	if (iret .ne. nf_noerr) goto 3601
466	if (ivtype .eq. 2) then
467	iret = nf_get_att_text( ncid, nf_global, attname, dumch1000 )
468	if (iret .ne. nf_noerr) goto 3601
469	i = max(1,min(1000,lenatt))
470
471	if (idiagaa .gt. 0) write(73,91010) &
472	iatt, attname(1:40), ivtype, lenatt, dumch1000(1:i)
473	else if (ivtype .eq. 4) then
474	iret = nf_get_att_int( ncid, nf_global, attname, idum )
475	if (iret .ne. nf_noerr) goto 3601
476	if (idiagaa .gt. 0) write(73,91020) &
477	iatt, attname(1:40), ivtype, lenatt, idum
478	else if ((ivtype .eq. 5) .and. (lenatt .eq. 1)) then
479	iret = nf_get_att_real( ncid, nf_global, attname, duma )
480	if (iret .ne. nf_noerr) goto 3601
481	if (idiagaa .gt. 0) write(73,91030) &
482	iatt, attname(1:40), ivtype, lenatt, duma
483	else if ((ivtype .eq. 5) .and. (lenatt .gt. 1)) then
484	allocate( duma_alloc(lenatt) )
485	iret = nf_get_att_real( ncid, nf_global, attname, duma_alloc )
486	if (iret .ne. nf_noerr) goto 3601
487	if (idiagaa .gt. 0) then
488	write(73,91010) iatt, attname(1:40), ivtype, lenatt
489	write(73,91040) (duma_alloc(i), i=1,lenatt)
490	end if
491
492	deallocate( duma_alloc )
493	else
494	if (idiagaa .gt. 0) write(*,'(i4,1x,a,2(1x,i6))') &
495	iatt, attname(1:40), ivtype, lenatt
496	goto 3401
497	endif
498
499	3401 continue
500
501	goto 3901
502
503	3601 write(*,9110) 'gridcheck-windfiled:error inquiring attribute',iatt,fnamenc
504	write(73,9110) 'gridcheck-windfiled:error inquiring attribute',iatt,fnamenc
505	stop
506
507	3901 continue
508
509
510	!c--------------------------------------------------
511	do i=1, nvars_tot
512	varid=i
513	iret = nf_inq_var ( ncid, varid, varname, xtype, &
514	ndims, dimids, natts)
515	if (iret .ne. nf_noerr) then
516	write (,) 'error in nf_inq_var'
517	stop
518	end if
519	!if (idiagaa.eq.1) then
520	!!write(75,*)'wind field inquire',ncid, varid, varname, xtype, &
521	!!ndims, dimids, natts,'fine'
522	!end if
523
524	xtypev(i)=xtype
525	varnamev(i)=varname
526	ndimsv(i)=ndims
527	do j=1,maxdim
528	dimidsm(i,j)=dimids(j)
529	end do
530	!if (idiagaa.eq.1) then
531	!!write(71,9031)i,varnamev(i),xtypev(i),ndims !WRITE diagnostic to file
532	!end if
533	end do
534
535
536
537
538	!c---------------- Inquiring about varnames and real fields ---------------------
539	nvar_exp_in_meteo_field=17
540	varnamev(1)='U' !horizontal wind
541	varnamev(2)='V' !horizontal wind
542	varnamev(3)='OMEGA' ! vertical wind pa/s
543	varnamev(4)='T' !temperature
544	varnamev(5)='Q' ! Specific humidity
545	varnamev(6)='PS' !surface pressure
546	varnamev(7)='CLDTOT' !total cloud cover
547	varnamev(8)='U10' !ten meters wind
548	varnamev(9)='TREFHT' !2m temperature
549	varnamev(10)='PRECL' !large scale precipitation
550	varnamev(11)='PRECC' !convective precipitation
551	varnamev(12)='SHFLX' !sensible heat fluxes
552	varnamev(13)='TAUX' !surface stress east-west
553	varnamev(14)='TAUY' !surface stress north-south
554	varnamev(15)='QREFHT' !speciifc humidity
555	varnamev(16)='SNOWHLND' !water equivakent snow depth !double check if must be also summed SNOWHICE
556	varnamev(17)='FSDS' !downwelling solar flux atsurface! double checck if this correct for stomata opeining parameterizations
557	!---------------- check and load file contents
558	do i=1,nvar_exp_in_meteo_field
559	varname=varnamev(i)
560	call check_variable(varname,fnamenc,maxdim,nf_float, &
561	id_var,ndims,id_dim,ierr,ncid)
562	if (ierr.lt.0) goto 100
563	lendim_exp=0
564	do j = 1, ndims
565	lendim_exp(j) = lendim(id_dim(j))
566	end do
567
568	call allocatedumarray(ndims,lendim_exp,maxdim,nf_float)
569	do j=1, ndims
570	istart(j) = 1 !ndims
571	icount(j) = lendim_exp(j)
572	end do
573	if (ndims.eq.1) then
574	iret = &
575	nf_get_vara_real( ncid, id_var, istart, icount, dumarray1D_real)
576	if (iret .ne. nf_noerr) then
577	write(*,9100) 'error inquiring var value ' // varname, fnamenc
578	ierr = -5
579	goto 100
580	end if
581	else if (ndims.eq.2) then
582	iret = &
583	nf_get_vara_real( ncid, id_var, istart, icount, dumarray2D_real)
584	if (iret .ne. nf_noerr) then
585	write(*,9100) 'error inquiring var value ' // varname, fnamenc
586	ierr = -5
587	goto 100
588	end if
589	else if (ndims.eq.3) then
590	iret = &
591	nf_get_vara_real( ncid, id_var, istart, icount, dumarray3D_real)
592	if (iret .ne. nf_noerr) then
593	write(*,9100) 'error inquiring var value ' // varname, fnamenc
594	ierr = -5
595	goto 100
596	end if
597	else if (ndims.eq.4) then
598	iret = &
599	nf_get_vara_real( ncid, id_var, istart, icount, dumarray4D_real)
600	if (iret .ne. nf_noerr) then
601	write(*,9100) 'error inquiring var value ' // varname, fnamenc
602	ierr = -5
603	goto 100
604	end if
605	end if
606	!!write(75,*)'gridcheck read variable', varname
607	!c------------- read number of vertical level
608	!c------------- seems that U, V, T, Q, OMEGA are co-located in CAM3.0/CAM4.0, see user's guide to thE NCAR CAM 3.0 page 38-
609	!c------------- so nwz and nuvz are the same here while etadot levels from surface to top are nwz+1=27 for CAM 4.+0!
610	if (varname.eq.'OMEGA') then
611	nwz=lendim_exp(3)
612	else if (varname.eq.'U') then
613	nuvz=lendim_exp(3)
614	end if
615
616	!c----------
617
618	end do
619
620
621	!200 continue
622
623	!C If desired, shift all grids by nxshift grid cells
624	!*************************************************
625
626	if (xglobal) then
627	call shift_field_0(oro,nxfield,ny)
628	call shift_field_0(lsm,nxfield,ny)
629	call shift_field_0(excessoro,nxfield,ny)
630	endif
631
632	! Output of grid info
633	!********************
634
635	write(,)
636	write(,)
637	write(*,'(a,2i7)') '# of vertical levels in NORESM data: ', &
638	nuvz+1,nwz
639	write(,)
640	write(*,'(a)') 'Mother domain:'
641	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
642	xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
643	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ', &
644	ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
645	write(,)
646
647	!C-------------- load vertical grid strucutre --------------
648	nvar_exp_in_meteo_field=6
649	varnamev(1)='P0' !omega levels
650	varnamev(2)='hyai' !etadot levels
651	varnamev(3)='hybi' !etadot levels
652	varnamev(4)='hyam' !omega levels
653	varnamev(5)='hybm' !omega levels
654	varnamev(6)='time' !omega levels
655	do i=1,nvar_exp_in_meteo_field
656	varname=varnamev(i)
657	call check_variable(varname,fnamenc,maxdim,nf_double, &
658	id_var,ndims,id_dim,ierr,ncid)
659	if (ierr.lt.0) goto 100
660	lendim_exp=0
661	do j = 1, ndims
662	lendim_exp(j) = lendim(id_dim(j))
663	end do
664	! if (ndims.eq.0) then !this alow teh use of a 1D array of 1 element for reading scalar variable (0-D)
665	! ndims=1
666	! lendim_exp(1)=1
667	! end if
668	call allocatedumarray(ndims,lendim_exp,maxdim,nf_double)
669	do j=1, ndims
670	istart(j) = 1 !ndims
671	icount(j) = lendim_exp(j)
672	end do
673
674	if (ndims.eq.0) then
675	iret = &
676	nf_get_vara_double( ncid, id_var, 1, 1, dumvar)
677	if (iret .ne. nf_noerr) then
678	write(*,9100) 'error inquiring var value ' // varname, fnamenc
679	ierr = -5
680	goto 100
681	end if
682
683	else if (ndims.eq.1) then
684	iret = &
685	nf_get_vara_double( ncid, id_var, istart, icount, dumarray1D)
686	if (iret .ne. nf_noerr) then
687	write(*,9100) 'error inquiring var value ' // varname, fnamenc
688	ierr = -5
689	goto 100
690	end if
691	end if
692
693	if (varname.eq.'hyai') then !etadot level 27 levels
694	netadot=nwz+1
695	do j=1, nwz+1
696	akm(nwz+2-j)=sngl(dumarray1D(j)*dumvar) ! dumvar contains P0
697	end do
698	else if (varname.eq.'hybi') then !etadot level 27 levels
699	do j=1, nwz+1
700	bkm(nwz+2-j)=sngl(dumarray1D(j))
701	end do
702	else if (varname.eq.'hyam') then !omega and U,V levels 26 levels, use 27 because added an "artificial" layer at the ground
703	do j=1, nuvz
704	akz(nuvz+2-j)=sngl(dumarray1D(j)*dumvar) !dumvar contains P0
705	end do
706	akz(1)=0.
707	else if (varname.eq.'hybm') then !omega and U,V levels 26 levels, use 27 because added an "artificial" layer at the ground
708	do j=1, nuvz
709	bkz(nuvz+2-j)=sngl(dumarray1D(j))
710	end do
711	bkz(1)=1.0
712	end if
713	end do
714
715	!************* Output of grid info!******************
716
717	write(,)
718	write(,)
719	write(*,'(a,2i7)') '# of vertical levels in NORESM data: ', &
720	nuvz+1,nwz
721	write(,)
722	write(*,'(a)') 'Mother domain:'
723	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
724	xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
725	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ', &
726	ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
727	write(,)
728	nuvz=nuvz+1
729	nz=nuvz
730
731	!c-- assign vertical discretization ------------------
732
733	if (nz.gt.nzmax) stop 'nzmax too small'
734	do i=1,nuvz
735	aknew(i)=akz(i)
736	bknew(i)=bkz(i)
737	end do
738
739	! THE DOUBLING IS NOT active in NORESM VERSION : Switch on following lines to use doubled vertical resolution
740	! this is not active in FLEXPART-NorESM !!!!
741	! *************************************************************
742	! nz=nuvz+nwz-1
743	! if (nz.gt.nzmax) stop 'nzmax too small'
744	! do 100 i=1,nwz
745	! aknew(2*(i-1)+1)=akm(i)
746	! 100 bknew(2*(i-1)+1)=bkm(i)
747	! do 110 i=2,nuvz
748	! aknew(2*(i-1))=akz(i)
749	! 10 bknew(2*(i-1))=bkz(i)
750	! !End doubled vertical resolution
751
752	!C-------------- load date --------------
753	nvar_exp_in_meteo_field=1
754	varnamev(1)='datesec' !date
755
756
757	do i=1,nvar_exp_in_meteo_field
758	varname=varnamev(i)
759	call check_variable(varname,fnamenc,maxdim,nf_int, &
760	id_var,ndims,id_dim,ierr,ncid)
761	if (ierr.lt.0) goto 100
762	lendim_exp=0
763	do j = 1, ndims
764	lendim_exp(j) = lendim(id_dim(j))
765	end do
766	call allocatedumarray(ndims,lendim_exp,maxdim,nf_int)
767	do j=1, ndims
768	istart(j) = 1 !ndims
769	icount(j) = lendim_exp(j)
770	end do
771
772	iret = &
773	nf_get_vara_int( ncid, id_var, istart, icount, dumarray1D_int)
774	if (iret .ne. nf_noerr) then
775	write(*,9100) 'error inquiring var value ' // varname, fnamenc
776	ierr = -5
777	goto 100
778	end if
779	end do
780
781	iret = nf_close( ncid )
782
783	!c----------- deallocation to free memory 12-2013 --------
784	if (allocated(dumarray4D)) deallocate(dumarray4D)
785	if (allocated(dumarray3D)) deallocate(dumarray3D)
786	if (allocated(dumarray2D)) deallocate(dumarray2D)
787	if (allocated(dumarray1D)) deallocate(dumarray1D)
788	if (allocated(dumarray4D_real)) deallocate(dumarray4D_real)
789	if (allocated(dumarray3D_real)) deallocate(dumarray3D_real)
790	if (allocated(dumarray2D_real)) deallocate(dumarray2D_real)
791	if (allocated(dumarray1D_real)) deallocate(dumarray1D_real)
792	if (allocated(dumarray4D_int)) deallocate(dumarray4D_int)
793	if (allocated(dumarray3D_int)) deallocate(dumarray3D_int)
794	if (allocated(dumarray2D_int)) deallocate(dumarray2D_int)
795	if (allocated(dumarray1D_int)) deallocate(dumarray1D_int)
796
797
798	! Determine the uppermost level for which the convection scheme shall be applied
799	! by assuming that there is no convection above 50 hPa (for standard SLP)
800	!*****************************************************************************
801
802	do i=1,nuvz-2
803	pint=akz(i)+bkz(i)*101325.
804	if (pint.lt.5000.) goto 96
805	end do
806	96 nconvlev=i
807	if (nconvlev.gt.nconvlevmax-1) then
808	nconvlev=nconvlevmax-1
809	write(,) 'Attention, convection only calculated up to ', &
810	akz(nconvlev)+bkz(nconvlev)*101325,' Pa'
811	endif
812
813	return
814
815	100 continue
816	print *,'error reading',varname,'erros code',ierr
817	stop
818
819	return
820	end

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

Download in other formats:

Original Format