38 SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg)
39 use machine,
only: kind_phys, kind_io4, kind_io8
44 integer,
intent(in) :: me, master, iflip, idate(4)
45 character(len=*),
intent(inout) :: errmsg
46 integer,
intent(inout) :: errflg
49 integer :: ncid, varid, ndims, hmx
50 integer :: i, j, k, n, ii, imon, klev
51 character :: fname*50, mn*2, vname*10
53 integer :: dimids(nf90_max_var_dims)
54 integer :: dimlen(nf90_max_var_dims)
61 if (me == master)
then
62 if ( iflip == 0 )
then
63 print *,
"GFS is top-down"
65 print *,
"GFS is bottom-up"
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'
86 fname=trim(
"aeroclim.m"//
'01'//
".nc")
88 if(.not.netcdf_check(nf90_open(fname , nf90_nowrite, ncid), &
89 errmsg, errflg,
'open '//trim(fname)))
then
93 vname = trim(specname(1))
95 if(.not.netcdf_check(nf90_inq_varid(ncid, vname, varid), &
96 errmsg, errflg,
'find id of '//trim(vname)//
' var'))
then
100 if(.not.netcdf_check(nf90_inquire_variable(ncid, varid, ndims=ndims, dimids=dimids), &
101 errmsg, errflg,
'inquire details about '//trim(vname)//
' var'))
then
105 if(.not.netcdf_check(nf90_inquire_dimension(ncid, dimids(i), len=dimlen(i)), &
106 errmsg, errflg,
'inquire details about dimension'))
then
117 print *,
'MERRA2 dim: ',dimlen(1:ndims)
122 if (.not.
allocated(aer_lat))
then
123 allocate(aer_lat(latsaer))
124 allocate(aer_lon(lonsaer))
129 if(.not.netcdf_check(nf90_inq_varid(ncid,
'lat', varid), &
130 errmsg, errflg,
'find id of lat var'))
then
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
139 if(.not.netcdf_check(nf90_inq_varid(ncid,
'lon', varid), &
140 errmsg, errflg,
'find id of lon var'))
then
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
148 if(.not.netcdf_check(nf90_close(ncid), errmsg, errflg,
'close '//trim(fname)))
then
154 SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg)
155 use machine,
only: kind_phys, kind_dbl_prec
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
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
171 integer,
allocatable :: invardims(:)
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))
187 CALL w3movdat(rinc,idat,jdat)
192 call w3doxdat(jdat,jdow,jdoy,jday)
193 rjday = jdoy + jdat(5) / 24.
194 IF (rjday < aer_time(1)) rjday = rjday+365.
198 if (rjday < aer_time(j))
then
204 if (n2 > 12) n2 = n2 -12
206 call read_netfaer(n1, iflip, 1, errmsg, errflg)
208 call read_netfaer(n2, iflip, 2, errmsg, errflg)
216 SUBROUTINE setindxaer(npts,dlat,jindx1,jindx2,ddy,dlon, &
217 iindx1,iindx2,ddx,me,master)
220 use aerclm_def,
only: aer_lat, jaero=>latsaer, &
221 aer_lon, iaero=>lonsaer
226 integer npts, jindx1(npts),jindx2(npts),iindx1(npts),iindx2(npts)
227 real(kind=kind_phys) dlat(npts),ddy(npts),dlon(npts),ddx(npts)
232 jindx2(j) = jaero + 1
234 if (dlat(j) < aer_lat(i))
then
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)))
251 iindx2(j) = iaero + 1
253 if (dlon(j) < aer_lon(i))
then
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)))
274 SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, &
275 ddy,iindx1,iindx2,ddx,lev,prsl,aerout, errmsg,errflg)
277 use machine,
only: kind_phys, kind_dbl_prec
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
290 integer jindx1(npts), jindx2(npts), iindx1(npts), iindx2(npts)
291 integer me,idate(4), master, nthrds
292 integer idat(8),jdat(8)
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
311 CALL w3movdat(rinc,idat,jdat)
316 call w3doxdat(jdat,jdow,jdoy,jday)
317 rjday = jdoy + jdat(5) / 24.
318 IF (rjday < aer_time(1)) rjday = rjday+365.
322 if (rjday < aer_time(j))
then
328 if (n2 > 12) n2 = n2 -12
332 if (me == master)
write(*,*)
"read in a new month MERRA2", n2
338 aerin(i,j,k,ii,1) = aerin(i,j,k,ii,2)
344 call read_netfaer(n2, iflip, 2, errmsg, errflg)
349 tx1 = (aer_time(n2) - rjday) / (aer_time(n2) - aer_time(n1))
351 if (n2 > 12) n2 = n2 -12
357 temiy(j) = temi*ddy(j)
358 temjx(j) = temj*ddx(j)
359 ddxy(j) = ddx(j)*ddy(j)
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))
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))
402 if(prsl(j,l) >= aerpres(j,1))
then
404 aerout(j,l,ii) = aerpm(j,1,ii)
406 else if(prsl(j,l) <= aerpres(j,levsaer))
then
408 aerout(j,l,ii) = aerpm(j,levsaer,ii)
412 IF(prsl(j,l) <= aerpres(j,k) .and. prsl(j,l)>aerpres(j,k+1))
then
414 i2 = min(k+1,levsaer)
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
422 aerout(j,l,ii) = aerpm(j,i1,ii)*tx1 + aerpm(j,i2,ii)*tx2
436 subroutine read_netfaer(nf, iflip,nt, errmsg,errflg)
437 use machine,
only: kind_phys, kind_io4, kind_io8
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
450 allocate (buff(lonsaer, latsaer, levsw))
451 allocate (pres_tmp(lonsaer, levsw))
452 allocate (buffx(lonsaer, latsaer, levsw, 1))
460 write(mn,
'(i2.2)') nf
461 fname=trim(
"aeroclim.m"//mn//
".nc")
463 if(.not.netcdf_check(nf90_open(fname , nf90_nowrite, ncid), &
464 errmsg, errflg,
'open '//trim(fname)))
then
470 if(.not.netcdf_check(nf90_inq_varid(ncid,
"DELP", varid), &
471 errmsg, errflg,
'find id of DELP var'))
then
474 if(.not.netcdf_check(nf90_get_var(ncid, varid, buff), &
475 errmsg, errflg,
'read DELP var'))
then
484 pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k)
490 if ( iflip == 0 )
then
493 klev = ( levsw - k ) + 1
496 aer_pres(i,j,k,nt) = 1.d0*pres_tmp(i,klev)
505 vname=trim(specname(ii))
507 if(.not.netcdf_check(nf90_inq_varid(ncid, vname, varid), &
508 errmsg, errflg,
'get id of '//trim(vname)//
' var'))
then
511 if(.not.netcdf_check(nf90_get_var(ncid, varid, buffx), &
512 errmsg, errflg,
'read '//trim(vname)//
' var'))
then
519 if ( iflip == 0 )
then
522 klev = ( levsw - k ) + 1
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
536 if(.not.netcdf_check(nf90_close(ncid), errmsg, errflg,
'close '//trim(fname)))
then
539 deallocate (buff, pres_tmp)