CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
aerinterp.F90
1
4
9
10 implicit none
11
12 private read_netfaer
13
14 public :: read_aerdata, setindxaer, aerinterpol,read_aerdataf
15
16contains
17
18 logical function netcdf_check(status, errmsg, errflg, why)
19 use netcdf
20 implicit none
21 character(len=*), intent(inout) :: errmsg
22 integer, intent(out) :: errflg
23 integer, intent(in) :: status
24 character(len=*), intent(in) :: why
25
26 netcdf_check = (status == nf90_noerr)
27
28 if(netcdf_check) then
29 errflg = 0
30 errmsg = ' '
31 else
32 errflg = 1
33 errmsg = trim(why) // ': ' // trim(nf90_strerror(status))
34 endif
35
36 END function netcdf_check
37
38 SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg)
39 use machine, only: kind_phys, kind_io4, kind_io8
40 use aerclm_def
41 use netcdf
42
43!--- in/out
44 integer, intent(in) :: me, master, iflip, idate(4)
45 character(len=*), intent(inout) :: errmsg
46 integer, intent(inout) :: errflg
47
48!--- locals
49 integer :: ncid, varid, ndims, hmx
50 integer :: i, j, k, n, ii, imon, klev
51 character :: fname*50, mn*2, vname*10
52 logical :: file_exist
53 integer :: dimids(nf90_max_var_dims)
54 integer :: dimlen(nf90_max_var_dims)
55
56 errflg = 0
57 errmsg = ' '
58
59!
60!! ===================================================================
61 if (me == master) then
62 if ( iflip == 0 ) then ! data from toa to sfc
63 print *, "GFS is top-down"
64 else
65 print *, "GFS is bottom-up"
66 endif
67 endif
68!
69!! ===================================================================
70!! check if all necessary files exist
71!! ===================================================================
72 do imon = 1, 12
73 write(mn,'(i2.2)') imon
74 fname=trim("aeroclim.m"//mn//".nc")
75 inquire (file = fname, exist = file_exist)
76 if (.not. file_exist) then
77 errmsg = 'Error in read_aerdata: file ' // trim(fname) // ' not found'
78 errflg = 1
79 return
80 endif
81 enddo
82!
83!! ===================================================================
84!! fetch dim spec and lat/lon from m01 data set
85!! ===================================================================
86 fname=trim("aeroclim.m"//'01'//".nc")
87 ncid = -1
88 if(.not.netcdf_check(nf90_open(fname , nf90_nowrite, ncid), &
89 errmsg, errflg, 'open '//trim(fname))) then
90 return
91 endif
92
93 vname = trim(specname(1))
94 varid = -1
95 if(.not.netcdf_check(nf90_inq_varid(ncid, vname, varid), &
96 errmsg, errflg, 'find id of '//trim(vname)//' var')) then
97 return
98 endif
99 ndims = 0
100 if(.not.netcdf_check(nf90_inquire_variable(ncid, varid, ndims=ndims, dimids=dimids), &
101 errmsg, errflg, 'inquire details about '//trim(vname)//' var')) then
102 return
103 endif
104 do i=1,ndims
105 if(.not.netcdf_check(nf90_inquire_dimension(ncid, dimids(i), len=dimlen(i)), &
106 errmsg, errflg, 'inquire details about dimension')) then
107 return
108 endif
109 enddo
110
111! specify latsaer, lonsaer, hmx
112 lonsaer = dimlen(1)
113 latsaer = dimlen(2)
114 levsw = dimlen(3)
115
116 if(me==master) then
117 print *, 'MERRA2 dim: ',dimlen(1:ndims)
118 endif
119
120! allocate arrays
121
122 if (.not. allocated(aer_lat)) then
123 allocate(aer_lat(latsaer))
124 allocate(aer_lon(lonsaer))
125 endif
126
127! construct lat/lon array
128 varid = -1
129 if(.not.netcdf_check(nf90_inq_varid(ncid, 'lat', varid), &
130 errmsg, errflg, 'find id of lat var')) then
131 return
132 endif
133 aer_lat = 0
134 if(.not.netcdf_check(nf90_get_var(ncid, varid, aer_lat, (/ 1, 1, 1 /), (/latsaer, 1, 1/)), &
135 errmsg, errflg, 'read lat var')) then
136 return
137 endif
138 varid = -1
139 if(.not.netcdf_check(nf90_inq_varid(ncid, 'lon', varid), &
140 errmsg, errflg, 'find id of lon var')) then
141 return
142 endif
143 aer_lon = 0
144 if(.not.netcdf_check(nf90_get_var(ncid, varid, aer_lon, (/ 1, 1, 1 /), (/lonsaer, 1, 1/)), &
145 errmsg, errflg, 'read lon var')) then
146 return
147 endif
148 if(.not.netcdf_check(nf90_close(ncid), errmsg, errflg, 'close '//trim(fname))) then
149 return
150 endif
151 END SUBROUTINE read_aerdata
152!
153!**********************************************************************
154 SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg)
155 use machine, only: kind_phys, kind_dbl_prec
156 use aerclm_def
157
158!--- in/out
159 integer, intent(in) :: me, master, iflip, idate(4)
160 character(len=*), intent(inout) :: errmsg
161 integer, intent(inout) :: errflg
162 real(kind=kind_phys), intent(in) :: fhour
163!--- locals
164 integer :: i, j, k, n, ii, imon, klev, n1, n2
165 logical :: file_exist
166 integer idat(8),jdat(8)
167 real(kind=kind_phys) rjday
168 real(kind=kind_dbl_prec) rinc(5)
169 integer jdow, jdoy, jday
170
171 integer, allocatable :: invardims(:)
172!
173 if (.not. allocated(aerin)) then
174 allocate(aerin(iamin:iamax,jamin:jamax,levsaer,ntrcaerm,timeaer))
175 allocate(aer_pres(iamin:iamax,jamin:jamax,levsaer,timeaer))
176 endif
177
178! allocate local working arrays
179!! found interpolation months
180 idat = 0
181 idat(1) = idate(4)
182 idat(2) = idate(2)
183 idat(3) = idate(3)
184 idat(5) = idate(1)
185 rinc = 0.
186 rinc(2) = fhour
187 CALL w3movdat(rinc,idat,jdat)
188!
189 jdow = 0
190 jdoy = 0
191 jday = 0
192 call w3doxdat(jdat,jdow,jdoy,jday)
193 rjday = jdoy + jdat(5) / 24.
194 IF (rjday < aer_time(1)) rjday = rjday+365.
195!
196 n2 = 13
197 do j=2, 12
198 if (rjday < aer_time(j)) then
199 n2 = j
200 exit
201 endif
202 enddo
203 n1 = n2 - 1
204 if (n2 > 12) n2 = n2 -12
205!! ===================================================================
206 call read_netfaer(n1, iflip, 1, errmsg, errflg)
207 if(errflg/=0) return
208 call read_netfaer(n2, iflip, 2, errmsg, errflg)
209 if(errflg/=0) return
210!! ===================================================================
211 n1sv=n1
212 n2sv=n2
213!---
214 END SUBROUTINE read_aerdataf
215!
216 SUBROUTINE setindxaer(npts,dlat,jindx1,jindx2,ddy,dlon, &
217 iindx1,iindx2,ddx,me,master)
218!
219 USE machine, ONLY: kind_phys
220 use aerclm_def, only: aer_lat, jaero=>latsaer, &
221 aer_lon, iaero=>lonsaer
222!
223 implicit none
224!
225 integer me, master
226 integer npts, jindx1(npts),jindx2(npts),iindx1(npts),iindx2(npts)
227 real(kind=kind_phys) dlat(npts),ddy(npts),dlon(npts),ddx(npts)
228!
229 integer i,j
230
231 DO j=1,npts
232 jindx2(j) = jaero + 1
233 do i=1,jaero
234 if (dlat(j) < aer_lat(i)) then
235 jindx2(j) = i
236 exit
237 endif
238 enddo
239 jindx1(j) = max(jindx2(j)-1,1)
240 jindx2(j) = min(jindx2(j),jaero)
241 if (jindx2(j) .ne. jindx1(j)) then
242 ddy(j) = (dlat(j) - aer_lat(jindx1(j))) &
243 / (aer_lat(jindx2(j)) - aer_lat(jindx1(j)))
244 else
245 ddy(j) = 1.0
246 endif
247
248 ENDDO
249
250 DO j=1,npts
251 iindx2(j) = iaero + 1
252 do i=1,iaero
253 if (dlon(j) < aer_lon(i)) then
254 iindx2(j) = i
255 exit
256 endif
257 enddo
258 iindx1(j) = max(iindx2(j)-1,1)
259 iindx2(j) = min(iindx2(j),iaero)
260 if (iindx2(j) .ne. iindx1(j)) then
261 ddx(j) = (dlon(j) - aer_lon(iindx1(j))) &
262 / (aer_lon(iindx2(j)) - aer_lon(iindx1(j)))
263 else
264 ddx(j) = 1.0
265 endif
266 ENDDO
267
268 RETURN
269 END SUBROUTINE setindxaer
270!
271!**********************************************************************
272!**********************************************************************
273!
274 SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, &
275 ddy,iindx1,iindx2,ddx,lev,prsl,aerout, errmsg,errflg)
276!
277 use machine, only: kind_phys, kind_dbl_prec
278 use aerclm_def
279
280 implicit none
281 integer, intent(inout) :: errflg
282 character(*), intent(inout) :: errmsg
283 integer, intent(in) :: iflip
284 integer i1,i2, iday,j,j1,j2,l,npts,nc,n1,n2,lev,k,i,ii, klev
285 real(kind=kind_phys) fhour,temj, tx1, tx2,temi, tem
286 real(kind=kind_phys), dimension(npts) :: temij,temiy,temjx,ddxy
287
288!
289
290 integer jindx1(npts), jindx2(npts), iindx1(npts), iindx2(npts)
291 integer me,idate(4), master, nthrds
292 integer idat(8),jdat(8)
293!
294 real(kind=kind_phys) ddy(npts), ddx(npts),ttt
295 real(kind=kind_phys) aerout(npts,lev,ntrcaer)
296 real(kind=kind_phys) aerpm(npts,levsaer,ntrcaer)
297 real(kind=kind_phys) prsl(npts,lev), aerpres(npts,levsaer)
298 real(kind=kind_phys) rjday
299 real(kind=kind_dbl_prec) rinc(5)
300 integer jdow, jdoy, jday
301!
302 errflg = 0
303 errmsg = ' '
304 idat = 0
305 idat(1) = idate(4)
306 idat(2) = idate(2)
307 idat(3) = idate(3)
308 idat(5) = idate(1)
309 rinc = 0.
310 rinc(2) = fhour
311 CALL w3movdat(rinc,idat,jdat)
312!
313 jdow = 0
314 jdoy = 0
315 jday = 0
316 call w3doxdat(jdat,jdow,jdoy,jday)
317 rjday = jdoy + jdat(5) / 24.
318 IF (rjday < aer_time(1)) rjday = rjday+365.
319!
320 n2 = 13
321 do j=2, 12
322 if (rjday < aer_time(j)) then
323 n2 = j
324 exit
325 endif
326 enddo
327 n1 = n2 - 1
328 if (n2 > 12) n2 = n2 -12
329! need to read a new month
330 if (n1.ne.n1sv) then
331#ifdef DEBUG
332 if (me == master) write(*,*)"read in a new month MERRA2", n2
333#endif
334 DO ii = 1, ntrcaerm
335 do j = jamin, jamax
336 do k = 1, levsaer
337 do i = iamin, iamax
338 aerin(i,j,k,ii,1) = aerin(i,j,k,ii,2)
339 enddo !i-loop (lon)
340 enddo !k-loop (lev)
341 enddo !j-loop (lat)
342 ENDDO ! ii-loop (ntracaerm)
343!! ===================================================================
344 call read_netfaer(n2, iflip, 2, errmsg, errflg)
345 n1sv=n1
346 n2sv=n2
347 end if
348!
349 tx1 = (aer_time(n2) - rjday) / (aer_time(n2) - aer_time(n1))
350 tx2 = 1.0 - tx1
351 if (n2 > 12) n2 = n2 -12
352
353 do j=1,npts
354 temj = 1.0 - ddy(j)
355 temi = 1.0 - ddx(j)
356 temij(j) = temi*temj
357 temiy(j) = temi*ddy(j)
358 temjx(j) = temj*ddx(j)
359 ddxy(j) = ddx(j)*ddy(j)
360 enddo
361
362#ifndef __GFORTRAN__
363!$OMP parallel num_threads(nthrds) default(none) &
364!$OMP shared(npts,ntrcaer,aerin,aer_pres,prsl) &
365!$OMP shared(ddx,ddy,jindx1,jindx2,iindx1,iindx2) &
366!$OMP shared(aerpm,aerpres,aerout,lev,nthrds) &
367!$OMP shared(temij,temiy,temjx,ddxy) &
368!$OMP private(l,j,k,ii,i1,i2,j1,j2,tem) &
369!$OMP copyin(tx1,tx2) firstprivate(tx1,tx2)
370
371!$OMP do
372#endif
373 DO l=1,levsaer
374 DO j=1,npts
375 j1 = jindx1(j)
376 j2 = jindx2(j)
377 i1 = iindx1(j)
378 i2 = iindx2(j)
379 DO ii=1,ntrcaer
380 aerpm(j,l,ii) = &
381 tx1*(temij(j)*aerin(i1,j1,l,ii,1)+ddxy(j)*aerin(i2,j2,l,ii,1) &
382 +temiy(j)*aerin(i1,j2,l,ii,1)+temjx(j)*aerin(i2,j1,l,ii,1))&
383 +tx2*(temij(j)*aerin(i1,j1,l,ii,2)+ddxy(j)*aerin(i2,j2,l,ii,2) &
384 +temiy(j)*aerin(i1,j2,l,ii,2)+temjx(j)*aerin(i2,j1,l,ii,2))
385 ENDDO
386
387 aerpres(j,l) = &
388 tx1*(temij(j)*aer_pres(i1,j1,l,1)+ddxy(j)*aer_pres(i2,j2,l,1) &
389 +temiy(j)*aer_pres(i1,j2,l,1)+temjx(j)*aer_pres(i2,j1,l,1))&
390 +tx2*(temij(j)*aer_pres(i1,j1,l,2)+ddxy(j)*aer_pres(i2,j2,l,2) &
391 +temiy(j)*aer_pres(i1,j2,l,2)+temjx(j)*aer_pres(i2,j1,l,2))
392 ENDDO
393 ENDDO
394#ifndef __GFORTRAN__
395!$OMP end do
396
397! don't flip, input is the same direction as GFS (bottom-up)
398!$OMP do
399#endif
400 DO j=1,npts
401 DO l=1,lev
402 if(prsl(j,l) >= aerpres(j,1)) then
403 DO ii=1, ntrcaer
404 aerout(j,l,ii) = aerpm(j,1,ii) !! sfc level
405 ENDDO
406 else if(prsl(j,l) <= aerpres(j,levsaer)) then
407 DO ii=1, ntrcaer
408 aerout(j,l,ii) = aerpm(j,levsaer,ii) !! toa top
409 ENDDO
410 else
411 DO k=1, levsaer-1 !! from sfc to toa
412 IF(prsl(j,l) <= aerpres(j,k) .and. prsl(j,l)>aerpres(j,k+1)) then
413 i1 = k
414 i2 = min(k+1,levsaer)
415 exit
416 ENDIF
417 ENDDO
418 tem = 1.0 / (aerpres(j,i1) - aerpres(j,i2))
419 tx1 = (prsl(j,l) - aerpres(j,i2)) * tem
420 tx2 = (aerpres(j,i1) - prsl(j,l)) * tem
421 DO ii = 1, ntrcaer
422 aerout(j,l,ii) = aerpm(j,i1,ii)*tx1 + aerpm(j,i2,ii)*tx2
423 ENDDO
424 endif
425 ENDDO !L-loop
426 ENDDO !J-loop
427#ifndef __GFORTRAN__
428!$OMP end do
429
430!$OMP end parallel
431#endif
432
433 RETURN
434 END SUBROUTINE aerinterpol
435
436 subroutine read_netfaer(nf, iflip,nt, errmsg,errflg)
437 use machine, only: kind_phys, kind_io4, kind_io8
438 use aerclm_def
439 use netcdf
440 integer, intent(in) :: iflip, nf, nt
441 integer, intent(inout) :: errflg
442 character(*), intent(inout) :: errmsg
443 integer :: ncid, varid, i,j,k,ii,klev
444 character :: fname*50, mn*2, vname*10
445 real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff
446 real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx
447 real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp
448
449!! ===================================================================
450 allocate (buff(lonsaer, latsaer, levsw))
451 allocate (pres_tmp(lonsaer, levsw))
452 allocate (buffx(lonsaer, latsaer, levsw, 1))
453
454 errflg = 0
455 errmsg = ' '
456 buff = 0
457 pres_tmp = 0
458 buffx = 0
459
460 write(mn,'(i2.2)') nf
461 fname=trim("aeroclim.m"//mn//".nc")
462 ncid = -1
463 if(.not.netcdf_check(nf90_open(fname , nf90_nowrite, ncid), &
464 errmsg, errflg, 'open '//trim(fname))) then
465 return
466 endif
467
468! ====> construct 3-d pressure array (Pa)
469 varid = -1
470 if(.not.netcdf_check(nf90_inq_varid(ncid, "DELP", varid), &
471 errmsg, errflg, 'find id of DELP var')) then
472 return
473 endif
474 if(.not.netcdf_check(nf90_get_var(ncid, varid, buff), &
475 errmsg, errflg, 'read DELP var')) then
476 return
477 endif
478
479 do j = jamin, jamax
480 do i = iamin, iamax
481! constract pres_tmp (top-down), note input is top-down
482 pres_tmp(i,1) = 0.
483 do k=2, levsw
484 pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k)
485 enddo !k-loop
486 enddo !i-loop (lon)
487
488! extract pres_tmp to fill aer_pres (in Pa)
489 do k = 1, levsaer
490 if ( iflip == 0 ) then ! data from toa to sfc
491 klev = k
492 else ! data from sfc to top
493 klev = ( levsw - k ) + 1
494 endif
495 do i = iamin, iamax
496 aer_pres(i,j,k,nt) = 1.d0*pres_tmp(i,klev)
497 enddo !i-loop (lon)
498 enddo !k-loop (lev)
499 enddo !j-loop (lat)
500
501! ====> construct 4-d aerosol array (kg/kg)
502! merra2 data is top down
503! for GFS, iflip 0: toa to sfc; 1: sfc to toa
504 DO ii = 1, ntrcaerm
505 vname=trim(specname(ii))
506 varid = -1
507 if(.not.netcdf_check(nf90_inq_varid(ncid, vname, varid), &
508 errmsg, errflg, 'get id of '//trim(vname)//' var')) then
509 return
510 endif
511 if(.not.netcdf_check(nf90_get_var(ncid, varid, buffx), &
512 errmsg, errflg, 'read '//trim(vname)//' var')) then
513 return
514 endif
515
516 do j = jamin, jamax
517 do k = 1, levsaer
518! input is from toa to sfc
519 if ( iflip == 0 ) then ! data from toa to sfc
520 klev = k
521 else ! data from sfc to top
522 klev = ( levsw - k ) + 1
523 endif
524 do i = iamin, iamax
525 aerin(i,j,k,ii,nt) = 1.d0*buffx(i,j,klev,1)
526 if(aerin(i,j,k,ii,nt) < 0 .or. aerin(i,j,k,ii,nt) > 1.) then
527 aerin(i,j,k,ii,nt) = 1.e-15
528 endif
529 enddo !i-loop (lon)
530 enddo !k-loop (lev)
531 enddo !j-loop (lat)
532
533 ENDDO ! ii-loop (ntracaerm)
534
535! close the file
536 if(.not.netcdf_check(nf90_close(ncid), errmsg, errflg, 'close '//trim(fname))) then
537 return
538 endif
539 deallocate (buff, pres_tmp)
540 deallocate (buffx)
541 END SUBROUTINE read_netfaer
542
543end module aerinterp
544
This module contain subroutines of reading and interpolating aerosol data for MG microphysics.
Definition aerinterp.F90:8