Line data Source code
1 : #define _FILE 'physics/cam1/boundarydata.F90 '
2 : module boundarydata
3 : use shr_kind_mod, only: r8 => shr_kind_r8
4 : use spmd_utils, only: masterproc
5 : use ppgrid, only: pcols, pver, begchunk, endchunk
6 : use physics_types, only: physics_state
7 : use cam_abortutils, only: endrun
8 : #if ( defined SPMD )
9 : use mpishorthand, only: mpicom, mpir8, mpiint
10 : #endif
11 : use netcdf
12 : use error_messages, only: handle_ncerr
13 : use cam_logfile, only: iulog
14 : implicit none
15 : private
16 : integer, parameter :: ptrtim=12, ptrlon=1
17 :
18 : type boundarydata_type
19 : integer :: ncid
20 : integer :: fieldcnt
21 : integer :: nm
22 : integer :: np
23 : integer :: latsiz
24 : integer :: levsiz
25 : integer :: ncolsiz
26 : integer :: timsiz
27 : integer :: vertextrap
28 : logical :: iszonal, isncol
29 : integer :: ndims
30 : integer :: thistimedim
31 : integer :: psid
32 : integer :: map(4)
33 : integer :: dimids(4)
34 : integer, pointer :: dataid(:) => null()
35 : integer, pointer :: columnmap(:) => null()
36 : integer, pointer :: start(:,:) => null()
37 : integer, pointer :: count(:,:) => null()
38 : real(r8), pointer :: lat(:) => null()
39 : real(r8), pointer :: zi(:) => null()
40 : real(r8), pointer :: pin(:) => null()
41 : real(r8), pointer :: cdates(:) => null()
42 : real(r8), pointer :: fields(:,:,:,:,:) => null()
43 : real(r8), pointer :: datainst(:,:,:,:) => null()
44 : real(r8), pointer :: hybi(:) => null()
45 : real(r8), pointer :: ps(:,:,:) => null()
46 : end type boundarydata_type
47 :
48 : public boundarydata_init
49 : public boundarydata_update
50 : public boundarydata_vert_interp
51 : public boundarydata_type
52 :
53 : contains
54 3072 : subroutine boundarydata_init(bndyfilename,phys_state,fieldnames,fieldcnt,bndydata,vertextrap)
55 : implicit none
56 : character(len=*),intent(in) :: bndyfilename
57 : type(physics_state), intent(in):: phys_state(begchunk:endchunk)
58 : integer,intent(in) :: fieldcnt
59 : character(len=*), intent(in) :: fieldnames(fieldcnt)
60 : type(boundarydata_type),intent(out) :: bndydata
61 : integer,intent(in), optional :: vertextrap ! if 0 set values outside output grid to 0
62 : ! if 1 set to boundary value
63 : ! if 2 set to cyclic boundaries
64 : ! if 3 leave on input data grid and extrapolate later
65 3072 : real(r8), pointer :: datain(:,:,:,:,:)
66 : integer :: lchnk
67 :
68 3072 : bndydata%fieldcnt=fieldcnt
69 3072 : if(present(vertextrap)) then
70 0 : bndydata%vertextrap=vertextrap
71 : else
72 3072 : bndydata%vertextrap=0
73 : end if
74 : nullify(bndydata%fields)
75 :
76 3072 : call boundarydata_read(phys_state,bndyfilename,fieldcnt,fieldnames,bndydata,datain)
77 :
78 3072 : if(bndydata%iszonal) then
79 3072 : call boundarydata_interpolate(phys_state,datain,bndydata)
80 :
81 : allocate(bndydata%datainst(size(bndydata%fields,1),size(bndydata%fields,2), &
82 18432 : begchunk:endchunk,bndydata%fieldcnt))
83 :
84 3072 : deallocate(datain)
85 : end if
86 3072 : end subroutine boundarydata_init
87 :
88 741888 : subroutine boundarydata_update(phys_state, bndydata, update_out)
89 : use interpolate_data,only : get_timeinterp_factors
90 : type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
91 : type(boundarydata_type), intent(inout) :: bndydata
92 : logical, intent(out), optional :: update_out
93 : real(r8) :: cdate
94 : integer :: nm, np, lchnk, j, k, fld, cols, cole, ncol, ndims
95 : real(r8) :: fact1, fact2
96 741888 : real(r8), allocatable :: datain(:,:,:,:,:)
97 : logical :: update
98 : integer :: latspan
99 : integer :: kmax
100 : integer :: count(4), start(4), ierr
101 :
102 :
103 741888 : call get_data_bounding_date_indices(bndydata%cdates,bndydata%nm,bndydata%np,cdate,update)
104 741888 : if(present(update_out)) update_out=update
105 741888 : nm= bndydata%nm
106 741888 : np= bndydata%np
107 :
108 741888 : call get_timeinterp_factors(.true., np, bndydata%cdates(nm), bndydata%cdates(np), &
109 1483776 : cdate, fact1, fact2, _FILE)
110 :
111 741888 : if(size(bndydata%fields,5).eq.2) then
112 741888 : nm=1
113 741888 : np=2
114 741888 : if(update) then ! we need to read in the next month and interpolate
115 0 : if(bndydata%isncol) then
116 0 : bndydata%fields(:,:,:,:,nm)=bndydata%fields(:,:,:,:,np)
117 0 : do lchnk=begchunk,endchunk
118 0 : ncol=phys_state(lchnk)%ncol
119 0 : cols=1
120 0 : cole=cols+bndydata%count(cols,lchnk)-1
121 0 : do while(cole<=ncol)
122 :
123 0 : if(bndydata%levsiz==1) then
124 0 : ndims=2
125 0 : start=(/bndydata%start(cols,lchnk),bndydata%np,-1,-1/)
126 0 : count=(/bndydata%count(cols,lchnk),1,-1,-1/)
127 : else
128 0 : ndims=3
129 0 : start=(/bndydata%start(cols,lchnk),bndydata%levsiz,bndydata%np,-1/)
130 0 : count=(/bndydata%count(cols,lchnk),1,1,-1/)
131 : end if
132 0 : do fld=1,bndydata%fieldcnt
133 0 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(fld) , &
134 0 : bndydata%fields(cols:cole,:,lchnk,fld,np), &
135 0 : start(1:ndims), count(1:ndims)),&
136 0 : _FILE,__LINE__)
137 :
138 : end do
139 0 : if(cols==ncol) exit
140 0 : cols=cols+bndydata%count(cols,lchnk)
141 0 : cole=cols+bndydata%count(cols,lchnk)-1
142 : end do
143 : end do
144 :
145 : else
146 0 : allocate(datain(ptrlon,bndydata%levsiz,bndydata%latsiz,1,bndydata%fieldcnt))
147 0 : if(masterproc) then
148 0 : count(1)=ptrlon
149 0 : count(2)=bndydata%levsiz
150 0 : count(3)=bndydata%latsiz
151 0 : count(4)=1
152 0 : start(1)=1
153 0 : start(2)=1
154 0 : start(3)=1
155 0 : start(4)=bndydata%np
156 0 : write(iulog,*) 'boundarydata reading data for month: ',bndydata%np
157 0 : do fld=1,bndydata%fieldcnt
158 0 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(fld), &
159 0 : datain(:,:,:,:,fld), start, count),_FILE,__LINE__)
160 : end do
161 : end if
162 : #ifdef SPMD
163 0 : call mpibcast (datain, bndydata%levsiz*bndydata%latsiz*1*bndydata%fieldcnt, mpir8, 0, mpicom, ierr)
164 : #endif
165 0 : bndydata%fields(:,:,:,:,nm) = bndydata%fields(:,:,:,:,np)
166 0 : call boundarydata_interpolate(phys_state,datain,bndydata)
167 0 : deallocate(datain)
168 : end if
169 : end if
170 : end if
171 741888 : kmax = size(bndydata%fields,2)
172 :
173 2596608 : do fld=1,bndydata%fieldcnt
174 10073448 : do lchnk=begchunk,endchunk
175 7476840 : if(bndydata%isncol) then
176 0 : latspan = phys_state(lchnk)%ncol
177 : else
178 7476840 : latspan=phys_state(lchnk)%ulatcnt
179 : end if
180 704677680 : do k=1,kmax
181 11409491205 : do j=1,latspan
182 0 : bndydata%datainst(j,k,lchnk,fld)=bndydata%fields(j,k,lchnk,fld,nm)*fact1 + &
183 11402014365 : bndydata%fields(j,k,lchnk,fld,np)*fact2
184 : end do
185 : end do
186 : end do
187 : end do
188 741888 : end subroutine boundarydata_update
189 :
190 :
191 3072 : subroutine boundarydata_read(phys_state,bndyfilename,fieldcnt,fieldnames,bndydata,datain)
192 : !-----------------------------------------------------------------------
193 : !
194 : ! Purpose:
195 : ! Do initial read of time-variant boundary dataset, containing
196 : ! 12 monthly fields as a function of latitude and pressure. Determine the two
197 : ! consecutive months between which the current date lies.
198 : !
199 : ! Method:
200 : !
201 : ! Author: NCAR CMS
202 : !-----------------------------------------------------------------------
203 : use ioFileMod, only : getfil
204 : use bnddyi_mod, only: bnddyi
205 :
206 : implicit none
207 : type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
208 : character(len=*),intent(in) :: bndyfilename
209 : integer,intent(in) :: fieldcnt
210 : character(len=*), intent(in) :: fieldnames(fieldcnt)
211 : type(boundarydata_type), intent(inout) :: bndydata
212 : real(r8), pointer :: datain(:,:,:,:,:) !
213 : !
214 : ! Local variables
215 : !
216 : integer :: londimid
217 : integer :: latdimid
218 : integer :: ncoldimid
219 : integer :: levdimid
220 : integer :: ilevdimid
221 : integer :: timdimid
222 : integer :: ndims
223 : integer :: dimlen
224 : integer :: ilevsiz
225 : integer :: ncolsiz
226 : character(len=nf90_max_name) :: dimname
227 :
228 :
229 : ! integer :: ncid ! netcdf id for file
230 : integer :: dateid ! netcdf id for date variable
231 : integer :: secid ! netcdf id for seconds variable
232 : integer :: lonid ! netcdf id for longitude variable
233 : integer :: ncolid ! netcdf id for longitude variable
234 : integer :: latid ! netcdf id for latitude variable
235 : integer :: levid ! netcdf id for level variable
236 : integer :: timid ! netcdf id for time variable
237 : integer :: hybid
238 :
239 : integer :: dataid ! netcdf id for data fields
240 :
241 : integer :: lonsiz ! size of longitude dimension on tracer dataset
242 : integer :: levsiz ! size of level dimension on tracer dataset
243 : integer :: latsiz ! size of latitude dimension on tracer dataset
244 :
245 : integer :: j,n,k,nt,id ! indices
246 : integer :: ki,ko,ji,jo ! indices
247 : integer :: date_tr(ptrtim), sec_tr(ptrtim)
248 :
249 : integer :: dimids(4), start(4), count(4)
250 3072 : integer, pointer :: columnmap(:)
251 : real(r8) :: calday ! current calendar day
252 3072 : real(r8), pointer :: pin(:)
253 3072 : real(r8), allocatable :: tmp_ps(:,:), tmp_fld(:,:,:)
254 : integer :: mincid,maxcid
255 : real(r8), allocatable, target:: lati(:)
256 : integer :: cols, cole
257 : integer :: ierr, dimcnt
258 : integer :: i, ncol, lchnk
259 : character(len=256) :: locfn ! netcdf local filename to open
260 :
261 : !
262 : !-----------------------------------------------------------------------
263 : !
264 : ! SPMD: Master reads dataset and does broadcast. All subsequent interpolation is
265 : ! done in every process. This is not required, one could remove this conditional
266 : ! and read the dataset independently on each task.
267 : !
268 3072 : if(masterproc) then
269 4 : write(iulog,*)'boundarydata_read: Reading from: ', trim(bndyfilename)
270 : #ifndef USE_MASTERPROC
271 : end if
272 : #endif
273 3072 : call getfil(bndyfilename, locfn)
274 : call handle_ncerr( nf90_open(locfn, 0, bndydata%ncid),&
275 3072 : _FILE,__LINE__)
276 :
277 : ! write(iulog,*)'boundarydata_read: NCOPN returns id ',bndydata%ncid,' for file ',trim(locfn)
278 : !
279 : !------------------------------------------------------------------------
280 : ! Read tracer data
281 : !------------------------------------------------------------------------
282 : !
283 : ! Get dimension info
284 : !
285 3072 : nullify(columnmap)
286 3072 : nullify(pin)
287 : call handle_ncerr( nf90_inquire(bndydata%ncid, bndydata%ndims, unlimiteddimid=timdimid), &
288 3072 : _FILE,__LINE__)
289 3072 : ncolsiz=-1
290 3072 : levsiz=-1
291 3072 : lonsiz=-1
292 3072 : latsiz=-1
293 3072 : bndydata%isncol = .false.
294 15360 : do i=1,bndydata%ndims
295 : call handle_ncerr( nf90_inquire_dimension(bndydata%ncid, i, dimname, dimlen),&
296 12288 : _FILE,__LINE__)
297 15360 : if (dimname(1:3).eq.'lat') then
298 3072 : latdimid=i
299 3072 : latsiz=dimlen
300 9216 : else if (dimname(1:3) .eq.'lon') then
301 3072 : londimid=i
302 3072 : lonsiz=dimlen
303 6144 : else if (dimname(1:4) .eq. 'ncol') then
304 0 : ncoldimid=i
305 0 : ncolsiz=dimlen
306 0 : bndydata%isncol=.true.
307 6144 : else if (dimname(1:3) .eq. 'lev') then
308 3072 : levdimid=i
309 3072 : levsiz=dimlen
310 3072 : else if (dimname(1:4) .eq. 'ilev') then
311 12288 : ilevdimid=i
312 12288 : ilevsiz=dimlen
313 3072 : else if (dimname(1:4) .eq. 'time') then
314 3072 : if(timdimid/=i) then
315 0 : timdimid=i
316 : end if
317 3072 : bndydata%timsiz=dimlen
318 : else
319 0 : write(iulog,*) 'Warning: do not know how to handle dimension ',&
320 0 : trim(dimname), ' in boundarydata.F90:313'
321 : end if
322 : end do
323 :
324 3072 : bndydata%iszonal = (latsiz>0 .and. lonsiz<=1)
325 3072 : if (bndydata%iszonal) then
326 9216 : allocate(bndydata%lat(latsiz))
327 : end if
328 3072 : if(bndydata%isncol) then
329 : ! allocate (columnmap(ncolsiz))
330 : ! call handle_ncerr( nf90_inq_varid(bndydata%ncid, 'ncol' , ncolid),&
331 : ! _FILE,__LINE__)
332 : ! call handle_ncerr( nf90_get_var(bndydata%ncid,ncolid,columnmap), &
333 : ! _FILE,__LINE__)
334 0 : if(levsiz>0) then
335 0 : allocate(bndydata%fields(pcols,levsiz,begchunk:endchunk,fieldcnt,2))
336 :
337 0 : ierr = nf90_inq_varid(bndydata%ncid, 'PS', bndydata%psid)
338 0 : if(ierr.eq.NF90_NOERR) then
339 0 : allocate(bndydata%ps(pcols,begchunk:endchunk,2))
340 0 : allocate(bndydata%hybi(levsiz+1))
341 : call handle_ncerr(nf90_inq_varid(bndydata%ncid,'hybi',hybid),&
342 0 : _FILE,__LINE__)
343 : call handle_ncerr( nf90_get_var(bndydata%ncid, hybid, bndydata%hybi ),&
344 0 : _FILE,__LINE__)
345 : else
346 0 : call endrun('BOUNDARYDATA_READ: Did not recognize a vertical coordinate variable')
347 : end if
348 : else
349 0 : levsiz=1
350 0 : allocate(bndydata%fields(pcols,1,begchunk:endchunk,fieldcnt,2))
351 : end if
352 : else
353 21504 : allocate(datain(lonsiz,levsiz,latsiz,2,fieldcnt))
354 : !
355 : ! Check dimension info
356 : !
357 3072 : if (lonsiz/=ptrlon) then
358 0 : call endrun ('BOUNDARYDATA_READ: longitude dependence not implemented')
359 : endif
360 :
361 3072 : if (bndydata%timsiz /= ptrtim) then
362 0 : write(iulog,*)'BOUNDARYDATA_READ: timsiz=',bndydata%timsiz,' must = ptrtim=',ptrtim
363 0 : call endrun
364 : end if
365 3072 : if( bndydata%vertextrap.lt.3) then
366 9216 : allocate(pin(levsiz))
367 : else
368 0 : allocate(bndydata%pin(levsiz))
369 0 : pin => bndydata%pin
370 : end if
371 9216 : allocate(bndydata%lat(latsiz))
372 :
373 :
374 18432 : allocate(datain(ptrlon,levsiz,latsiz,2,fieldcnt))
375 :
376 : call handle_ncerr( nf90_inq_varid(bndydata%ncid, 'lat' , latid),&
377 3072 : _FILE,__LINE__)
378 : end if
379 : !
380 : ! Determine necessary dimension and variable id's
381 : !
382 9216 : allocate(bndydata%cdates(bndydata%timsiz))
383 :
384 : call handle_ncerr(nf90_inq_varid(bndydata%ncid, 'date' , dateid), &
385 3072 : _FILE,__LINE__)
386 : call handle_ncerr( nf90_get_var(bndydata%ncid, dateid, date_tr),&
387 3072 : _FILE,__LINE__)
388 3072 : ierr = nf90_inq_varid(bndydata%ncid, 'datesec', secid)
389 3072 : if(ierr==NF90_NOERR) then
390 : call handle_ncerr( nf90_get_var(bndydata%ncid, secid , sec_tr),&
391 3072 : _FILE,__LINE__)
392 : else
393 0 : sec_tr=0
394 : end if
395 :
396 3072 : if (mod(date_tr(1),10000)/100 /= 1) then
397 0 : call endrun ('(boundarydata_read): error when cycling data: 1st month must be 1')
398 : end if
399 3072 : if (mod(date_tr(ptrtim),10000)/100 /= 12) then
400 0 : call endrun ('(boundarydata_read): error when cycling data: last month must be 12')
401 : end if
402 : !
403 : ! return the calander dates of the file data
404 : !
405 39936 : do n=1,ptrtim
406 39936 : call bnddyi(date_tr(n), sec_tr(n), bndydata%cdates(n))
407 : end do
408 : ! else
409 : ! call handle_ncerr( nf90_inq_varid(bndydata%ncid, 'time', dateid),&
410 : ! _FILE,__LINE__)
411 : !
412 : ! call handle_ncerr( nf90_get_var(bndydata%ncid, dateid, bndydata%cdates),&
413 : ! _FILE,__LINE__)
414 : !
415 : ! end if
416 : #ifdef USE_MASTERPROC
417 : else
418 : allocate(bndydata%cdates(ptrtim))
419 : end if
420 : #ifdef SPMD
421 : call mpibcast (bndydata%cdates, ptrtim, mpir8, 0, mpicom, ierr)
422 : #endif
423 : #endif
424 3072 : bndydata%nm=12
425 3072 : bndydata%np=1
426 3072 : call get_data_bounding_date_indices(bndydata%cdates,bndydata%nm,bndydata%np)
427 : #ifdef USE_MASTERPROC
428 : if(masterproc) then
429 : #endif
430 : !
431 : ! Obtain entire date and sec variables. Assume that will always
432 : ! cycle over 12 month data.
433 : !
434 : !
435 : ! Obtain input data latitude and level arrays.
436 : !
437 3072 : if(bndydata%iszonal) then
438 : call handle_ncerr( nf90_get_var(bndydata%ncid, latid, bndydata%lat),&
439 3072 : _FILE,__LINE__)
440 3072 : ierr = nf90_inq_varid(bndydata%ncid, 'lev' , levid)
441 : call handle_ncerr( nf90_get_var(bndydata%ncid, levid, pin ),&
442 3072 : _FILE,__LINE__)
443 : end if
444 :
445 9216 : allocate(bndydata%dataid(fieldcnt))
446 3072 : if(masterproc) then
447 4 : write(iulog,*) 'boundarydata reading data for months: ',bndydata%nm,bndydata%np
448 : end if
449 10752 : do i=1,fieldcnt
450 7680 : call handle_ncerr( nf90_inq_varid(bndydata%ncid, fieldnames(i) , bndydata%dataid(i)),&
451 18432 : _FILE,__LINE__)
452 : end do
453 3072 : if(bndydata%isncol) then
454 : allocate(bndydata%start(pcols,begchunk:endchunk), &
455 0 : bndydata%count(pcols,begchunk:endchunk))
456 :
457 : !
458 : ! For i/o efficiency we read in a block of data which includes the data needed on this
459 : ! processor but which may in fact include data not needed here. physics cids are just the
460 : ! offset into the file.
461 : !
462 :
463 0 : bndydata%start=-1
464 0 : bndydata%count=1
465 0 : mincid=2147483647
466 0 : maxcid=-1
467 0 : do lchnk=begchunk,endchunk
468 0 : ncol=phys_state(lchnk)%ncol
469 0 : i=minval(phys_state(lchnk)%cid(1:ncol))
470 0 : if(i < mincid) mincid = i
471 0 : i=maxval(phys_state(lchnk)%cid(1:ncol))
472 0 : if(i > maxcid) maxcid = i
473 : end do
474 :
475 0 : allocate(tmp_ps(mincid:maxcid,2))
476 0 : start=(/mincid,bndydata%nm,1,-1/)
477 0 : if(bndydata%np>bndydata%nm) then
478 0 : count=(/maxcid-mincid+1,2,-1,-1/)
479 : else
480 0 : count=(/maxcid-mincid+1,1,-1,-1/)
481 : end if
482 0 : if(associated(bndydata%ps) ) then
483 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%psid , &
484 0 : tmp_ps(:,1:count(2)), start(1:2), &
485 : count(1:2)),&
486 0 : _FILE,__LINE__)
487 0 : if(bndydata%np<bndydata%nm) then
488 0 : start(2)=bndydata%np
489 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%psid , &
490 : tmp_ps(:,2:2), start(1:2), &
491 : count(1:2)),&
492 0 : _FILE,__LINE__)
493 : end if
494 :
495 0 : do lchnk=begchunk,endchunk
496 0 : do n=1,phys_state(lchnk)%ncol
497 0 : bndydata%ps(n ,lchnk,:) = tmp_ps(phys_state(lchnk)%cid(n),:)
498 : end do
499 : end do
500 0 : deallocate(tmp_ps)
501 :
502 : end if
503 :
504 0 : if(levsiz>1) then
505 : dimcnt=3
506 : else
507 0 : dimcnt=2
508 : end if
509 0 : start(2)=1
510 0 : count(2)=levsiz
511 :
512 0 : if(bndydata%np>bndydata%nm) then
513 0 : count(dimcnt)=2
514 : else
515 0 : count(dimcnt)=1
516 : end if
517 0 : start(dimcnt)=bndydata%nm
518 :
519 0 : allocate(tmp_fld(mincid:maxcid,count(2),2))
520 :
521 0 : do i=1,fieldcnt
522 0 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i) , &
523 0 : tmp_fld(:,:,1:count(dimcnt)), &
524 : start(1:dimcnt), count(1:dimcnt)),&
525 0 : _FILE,__LINE__)
526 :
527 0 : do lchnk=begchunk,endchunk
528 0 : do n=1,phys_state(lchnk)%ncol
529 0 : bndydata%fields(n,:,lchnk,i,1:count(dimcnt)) = tmp_fld(phys_state(lchnk)%cid(n),:,:)
530 : end do
531 : end do
532 : end do
533 0 : if(bndydata%np<bndydata%nm) then
534 0 : start(dimcnt)=bndydata%np
535 0 : do i=1,fieldcnt
536 0 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i) , &
537 : tmp_fld(:,:,2:2), &
538 : start(1:dimcnt), count(1:dimcnt)),&
539 0 : _FILE,__LINE__)
540 0 : do lchnk=begchunk,endchunk
541 0 : do n=1,phys_state(lchnk)%ncol
542 0 : bndydata%fields(n,:,lchnk,i,1:count(dimcnt)) = tmp_fld(phys_state(lchnk)%cid(n),:,:)
543 : end do
544 : end do
545 : end do
546 : end if
547 0 : deallocate(tmp_fld)
548 : ! deallocate(columnmap)
549 : else
550 : !
551 : ! get the dimension orientation info from the first variable assume but verify that
552 : ! all variables requested have the same orientation
553 : !
554 3072 : allocate(bndydata%start(4,1),bndydata%count(4,1))
555 0 : call handle_ncerr( nf90_inquire_variable(bndydata%ncid,bndydata%dataid(1), &
556 3072 : ndims=bndydata%ndims,dimids=bndydata%dimids),_FILE,__LINE__)
557 :
558 18432 : bndydata%start=1
559 15360 : do id=1,bndydata%ndims
560 15360 : if(bndydata%dimids(id).eq.londimid) then
561 3072 : bndydata%map(id)=1
562 3072 : bndydata%count(id,1)=lonsiz
563 9216 : else if(bndydata%dimids(id).eq.levdimid) then
564 3072 : bndydata%map(id)=lonsiz
565 3072 : bndydata%count(id,1)=levsiz
566 6144 : else if(bndydata%dimids(id).eq.latdimid) then
567 3072 : bndydata%map(id)=lonsiz
568 6144 : if(any(bndydata%dimids.eq.levdimid)) bndydata%map(id)=bndydata%map(id)*levsiz
569 3072 : bndydata%count(id,1)=latsiz
570 3072 : else if(bndydata%dimids(id).eq.timdimid) then
571 3072 : bndydata%map(id)=lonsiz*latsiz
572 6144 : if(any(bndydata%dimids.eq.levdimid)) bndydata%map(id)=bndydata%map(id)*levsiz
573 3072 : bndydata%count(id,1)=2
574 3072 : bndydata%start(id,1)=bndydata%nm
575 3072 : bndydata%thistimedim=id
576 : else
577 0 : write(iulog,*) __LINE__,fieldnames(i),id,bndydata%dimids(id),londimid, &
578 0 : levdimid,latdimid,timdimid
579 0 : call endrun(_FILE)
580 : end if
581 : end do
582 :
583 10752 : do i=1,fieldcnt
584 7680 : call handle_ncerr( nf90_inq_varid(bndydata%ncid, fieldnames(i), bndydata%dataid(i)),&
585 15360 : _FILE,__LINE__)
586 :
587 0 : call handle_ncerr( nf90_inquire_variable(bndydata%ncid,bndydata%dataid(i), &
588 7680 : ndims=ndims,dimids=dimids),_FILE,__LINE__)
589 : if(ndims/=bndydata%ndims .or. dimids(1)/=bndydata%dimids(1).or.&
590 7680 : dimids(2)/=bndydata%dimids(2) .or. dimids(3)/=bndydata%dimids(3)) then
591 0 : call endrun('BOUNDARYDATA_READ: Variable dims or order does not match')
592 : end if
593 :
594 10752 : if(bndydata%np .gt. bndydata%nm) then
595 0 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i), &
596 : datain(:,:,:,:,i),bndydata%start(:,1), bndydata%count(:,1), &
597 0 : map=bndydata%map),_FILE,__LINE__)
598 : else
599 7680 : bndydata%count(bndydata%thistimedim,1)=1
600 7680 : bndydata%start(bndydata%thistimedim,1)=bndydata%nm
601 0 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i), &
602 0 : datain(:,:,:,1:1,i), bndydata%start(:,1), bndydata%count(:,1), &
603 7680 : map=bndydata%map), _FILE,__LINE__)
604 :
605 7680 : bndydata%start(bndydata%thistimedim,1)=bndydata%np
606 0 : call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i), &
607 0 : datain(:,:,:,2:2,i), bndydata%start(:,1), bndydata%count(:,1), &
608 7680 : map=bndydata%map), _FILE,__LINE__)
609 :
610 : endif
611 : end do
612 :
613 : end if
614 : #ifdef USE_MASTERPROC
615 : end if
616 : #ifdef SPMD
617 : call mpibcast (levsiz, 1, mpiint, 0, mpicom, ierr)
618 : call mpibcast (latsiz, 1, mpiint, 0, mpicom, ierr)
619 : #endif
620 : #endif
621 3072 : bndydata%latsiz=latsiz
622 3072 : bndydata%levsiz=levsiz
623 : #ifdef USE_MASTERPROC
624 : if(.not. masterproc) then
625 : if( bndydata%vertextrap.lt.3) then
626 : allocate(pin(levsiz))
627 : else
628 : allocate(bndydata%pin(levsiz))
629 : pin => bndydata%pin
630 : end if
631 : allocate(bndydata%lat(latsiz))
632 : allocate(datain(ptrlon,levsiz,latsiz,2,fieldcnt))
633 : endif
634 : #ifdef SPMD
635 : call mpibcast (bndydata%lat, latsiz, mpir8, 0, mpicom, ierr)
636 : call mpibcast (pin, levsiz, mpir8, 0, mpicom, ierr)
637 : call mpibcast (datain, levsiz*latsiz*2*fieldcnt, mpir8, 0, mpicom, ierr)
638 :
639 : #endif
640 : #endif
641 : ! Convert input pressure from millibars to pascals.
642 3072 : if(associated(pin)) then
643 175104 : pin=pin*100._r8
644 3072 : if(bndydata%vertextrap.lt.3) then
645 9216 : allocate(bndydata%zi(levsiz))
646 : !
647 : !
648 : ! Convert input pressure levels to height (m).
649 : !
650 175104 : do k=1,levsiz
651 175104 : bndydata%zi(k) = 7.0e3_r8 * log (1.0e5_r8 / pin(k))
652 : end do
653 3072 : deallocate(pin)
654 : end if
655 : end if
656 6144 : end subroutine boundarydata_read
657 :
658 3072 : subroutine boundarydata_interpolate(phys_state, datain, bndydata)
659 : use ref_pres, only : pref_mid
660 : use interpolate_data,only : interp_type, lininterp_init, &
661 : lininterp_finish, lininterp
662 : use physconst, only: pi
663 :
664 : type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
665 : real(r8),intent(in) :: datain(:,:,:,:,:)
666 : type(boundarydata_type), intent(inout) :: bndydata
667 : type(interp_type) :: interp_wgts, lev_wgts
668 :
669 : integer :: k, lchnk, nt, j, fcnt
670 : real(r8) :: zo(pver)
671 : real(r8) :: lato(pcols)
672 : integer :: ulatcnt
673 : integer :: maxlatcnt
674 : integer :: timesize, tvalout
675 :
676 : !------------------------------------------------------------------------
677 : ! Interpolate tracer data to model grid
678 : !------------------------------------------------------------------------
679 : !
680 : ! Loop over all input times.
681 : !
682 :
683 3072 : timesize=2
684 :
685 3072 : maxlatcnt=1
686 15456 : do lchnk=begchunk,endchunk
687 15456 : maxlatcnt=max(maxlatcnt,phys_state(lchnk)%ulatcnt)
688 : end do
689 3072 : if(bndydata%vertextrap.lt.3) then
690 : !
691 : ! Convert approximate cam pressure levels to height (m).
692 : !
693 288768 : do k=1,pver
694 288768 : zo (k) = 7.0e3_r8 * log (1.0e5_r8 / pref_mid(k))
695 : end do
696 :
697 3072 : call lininterp_init(bndydata%zi,size(bndydata%zi),zo,pver,bndydata%vertextrap,lev_wgts)
698 3072 : if(.not. associated(bndydata%fields)) then
699 18432 : allocate(bndydata%fields(maxlatcnt,pver,begchunk:endchunk,bndydata%fieldcnt,timesize))
700 97890876 : bndydata%fields=0_r8
701 : end if
702 : else
703 0 : if(.not. associated(bndydata%fields)) then
704 0 : allocate(bndydata%fields(maxlatcnt,bndydata%levsiz,begchunk:endchunk,bndydata%fieldcnt,timesize))
705 0 : bndydata%fields=0_r8
706 : end if
707 : endif
708 15456 : do lchnk=begchunk,endchunk
709 12384 : ulatcnt=phys_state(lchnk)%ulatcnt
710 :
711 : !
712 : ! Convert cam model latitudes to degrees.
713 : ! Input model latitudes already in degrees.
714 : !
715 203068 : do j=1,ulatcnt
716 203068 : lato(j) = phys_state(lchnk)%ulat(j)*180._r8/pi
717 : end do
718 :
719 12384 : call lininterp_init(bndydata%lat,size(bndydata%lat),lato(1:ulatcnt),ulatcnt,1,interp_wgts)
720 12384 : timesize = size(datain,4)
721 43344 : do fcnt=1,bndydata%fieldcnt
722 105264 : do nt = 1,timesize
723 61920 : if(timesize.gt.1) then
724 61920 : tvalout=nt
725 : else
726 : tvalout=2
727 : end if
728 92880 : if(bndydata%vertextrap.lt.3) then
729 : call lininterp(transpose(datain(1,:,:,nt,fcnt)),bndydata%latsiz,bndydata%levsiz, &
730 26916600 : bndydata%fields(1:ulatcnt,:,lchnk,fcnt,tvalout), ulatcnt, pver, interp_wgts, lev_wgts)
731 : else
732 0 : do k=1,bndydata%levsiz
733 0 : call lininterp(datain(1,k,:,nt,fcnt),bndydata%latsiz, &
734 0 : bndydata%fields(1:ulatcnt,k,lchnk,fcnt,tvalout), ulatcnt, interp_wgts)
735 : end do
736 : end if
737 : end do
738 : end do ! end loop over time samples
739 15456 : call lininterp_finish(interp_wgts)
740 : end do
741 3072 : if(bndydata%vertextrap.lt.3) &
742 3072 : call lininterp_finish(lev_wgts)
743 :
744 3072 : return
745 3072 : end subroutine boundarydata_interpolate
746 :
747 744960 : subroutine get_data_bounding_date_indices(cdates,nm,np, cdayout, update)
748 3072 : use time_manager, only: get_curr_date, get_perp_date, get_curr_calday, &
749 : is_perpetual
750 : real(r8), intent(in) :: cdates(ptrtim)
751 : real(r8), intent(out),optional :: cdayout
752 : logical, intent(out) ,optional :: update
753 : integer, intent(inout) :: nm, np
754 : integer :: n, np1
755 : real(r8) :: calday
756 : integer :: yr, mon, day ! components of a date
757 : integer :: ncdate ! current date in integer format [yyyymmdd]
758 : integer :: ncsec ! current time of day [seconds]
759 :
760 744960 : calday = get_curr_calday()
761 744960 : if(present(cdayout)) cdayout=calday
762 744960 : if(present(update)) update=.false. ! initialize output variable
763 :
764 : if(min(nm,np) .ge. 1 .and. max(nm,np) .le. 12 .and. &
765 744960 : calday>cdates(nm) .and. calday<=cdates(np)) return
766 744960 : if((nm==12 .and. np==1) .and. (calday <= cdates(np) .or. &
767 : calday > cdates(nm))) return
768 :
769 0 : if(present(update)) update=.true.
770 :
771 0 : if(calday <= cdates(1) .or. calday > cdates(12)) then
772 0 : nm=12
773 0 : np=1
774 : else
775 0 : nm=0
776 0 : do n=1,ptrtim-1
777 0 : if(calday>cdates(n) .and. calday<=cdates(n+1)) then
778 0 : nm=n
779 0 : np=n+1
780 : end if
781 : end do
782 0 : if(nm .eq. 0) then
783 0 : if ( is_perpetual() ) then
784 0 : call get_perp_date(yr, mon, day, ncsec)
785 : else
786 0 : call get_curr_date(yr, mon, day, ncsec)
787 : end if
788 0 : ncdate = yr*10000 + mon*100 + day
789 :
790 0 : write(iulog,*)'model date:', ncdate, ncsec,'boundary data dates:', cdates
791 0 : call endrun('get_data_bounding_date_indices: Failed to find dates bracketing dates')
792 : end if
793 : end if
794 :
795 744960 : end subroutine get_data_bounding_date_indices
796 :
797 :
798 : !================================================================================================
799 0 : subroutine boundarydata_vert_interp(lchnk, ncol, levsiz, fldcnt, pin, pmid, datain, dataout)
800 : !-----------------------------------------------------------------------
801 : !
802 : ! Purpose: Interpolate ozone from current time-interpolated values to model levels
803 : !
804 : ! Method: Use pressure values to determine interpolation levels
805 : !
806 : ! Author: Bruce Briegleb
807 : !
808 : !--------------------------------------------------------------------------
809 : implicit none
810 : ! Arguments
811 : !
812 : integer, intent(in) :: lchnk ! chunk identifier
813 : integer, intent(in) :: ncol ! number of atmospheric columns
814 : integer, intent(in) :: levsiz
815 : integer, intent(in) :: fldcnt
816 : real(r8), intent(in) :: pin(levsiz)
817 : real(r8), intent(in) :: pmid(pcols,pver) ! level pressures (mks)
818 : real(r8), intent(in) :: datain(pcols,levsiz,fldcnt)
819 : real(r8), intent(out) :: dataout(pcols,pver,fldcnt) ! ozone mass mixing ratio
820 : !
821 : ! local storage
822 : !
823 :
824 : integer :: i ! longitude index
825 : integer :: k, kk, kkstart ! level indices
826 : integer :: kupper(pcols) ! Level indices for interpolation
827 : integer :: kount ! Counter
828 : integer :: fld
829 : real(r8) dpu ! upper level pressure difference
830 : real(r8) dpl ! lower level pressure difference
831 : !--------------------------------------------------------------------------
832 : !
833 : ! Initialize index array
834 : !
835 0 : do i=1,ncol
836 0 : kupper(i) = 1
837 : end do
838 :
839 0 : do k=1,pver
840 : !
841 : ! Top level we need to start looking is the top level for the previous k
842 : ! for all longitude points
843 : !
844 : kkstart = levsiz
845 0 : do i=1,ncol
846 0 : kkstart = min0(kkstart,kupper(i))
847 : end do
848 : kount = 0
849 : !
850 : ! Store level indices for interpolation
851 : !
852 0 : do kk=kkstart,levsiz-1
853 0 : do i=1,ncol
854 0 : if (pin(kk).lt.pmid(i,k) .and. pmid(i,k).le.pin(kk+1)) then
855 0 : kupper(i) = kk
856 0 : kount = kount + 1
857 : end if
858 : end do
859 : !
860 : ! If all indices for this level have been found, do the interpolation and
861 : ! go to the next level
862 : !
863 0 : if (kount.eq.ncol) then
864 0 : do fld=1,fldcnt
865 0 : do i=1,ncol
866 0 : dpu = pmid(i,k) - pin(kupper(i))
867 0 : dpl = pin(kupper(i)+1) - pmid(i,k)
868 0 : dataout(i,k,fld) = (datain(i,kupper(i),fld )*dpl + &
869 0 : datain(i,kupper(i)+1,fld)*dpu)/(dpl + dpu)
870 :
871 : end do
872 : end do
873 : goto 35
874 : end if
875 : end do
876 : !
877 : ! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and
878 : ! must extrapolate from the bottom or top data level for at least some
879 : ! of the longitude points.
880 : !
881 0 : do fld=1,fldcnt
882 0 : do i=1,ncol
883 0 : if (pmid(i,k) .lt. pin(1)) then
884 0 : dataout(i,k,fld) = datain(i,1,fld)*pmid(i,k)/pin(1)
885 0 : else if (pmid(i,k) .gt. pin(levsiz)) then
886 0 : dataout(i,k,fld) = datain(i,levsiz,fld)
887 : else
888 0 : dpu = pmid(i,k) - pin(kupper(i))
889 0 : dpl = pin(kupper(i)+1) - pmid(i,k)
890 0 : dataout(i,k,fld) = (datain(i,kupper(i),fld )*dpl + &
891 0 : datain(i,kupper(i)+1,fld)*dpu)/(dpl + dpu)
892 : end if
893 : end do
894 : end do
895 0 : if (kount.gt.ncol) then
896 0 : call endrun ('ozone_data_vert_interp: Bad ozone data: non-monotonicity suspected')
897 : end if
898 0 : 35 continue
899 : end do
900 744960 : end subroutine boundarydata_vert_interp
901 : #if 0
902 : subroutine ncol_read_bracket(cid,columnmap,start,count,ncol)
903 : integer, intent(in) :: cid(:), columnmap(:), ncol
904 : integer, intent(out) :: start(:), count(:)
905 :
906 : integer :: i, j, tcol
907 :
908 : tcol = size(columnmap)
909 : count=1
910 : do i=1,ncol
911 : #if 1
912 : start(i)=cid(i)
913 : count(i)=1
914 : #else
915 : do j=1,tcol
916 : if(columnmap(j).eq.cid(i)) then
917 : start(i)=j
918 : count(i)=1
919 : exit
920 : end if
921 : end do
922 : #endif
923 : end do
924 : do i=1,ncol-1
925 : do j=1,ncol-i
926 : if(columnmap(start(i+j)).eq.columnmap(start(i)+j)) then
927 : count(i)=count(i)+1
928 : else
929 : exit
930 : end if
931 : end do
932 : end do
933 : write(iulog,*) __LINE__,cid(1),cid(ncol),minval(cid(1:ncol)),maxval(cid(1:ncol))
934 :
935 : end subroutine ncol_read_bracket
936 : #endif
937 0 : end module boundarydata
|
