94 function load_o3prog(this, file, fileID)
result (err_message)
96 integer,
intent(in) :: fileid
97 character(len=*),
intent(in) :: file
98 character(len=128) :: err_message
99 integer :: i1, i2, i3, ierr
100 real(kind=4), dimension(:),
allocatable :: lat4, pres4, time4, tempin
101 real(kind=4) :: blatc4
107 open(unit=fileid,file=trim(file), form=
'unformatted', convert=
'big_endian', iostat=ierr, iomsg=err_message)
108 if (ierr /= 0 )
return
109 read (fileid, iostat=ierr, iomsg=err_message) this%ncf, this%nlat, this%nlev, this%ntime
110 if (ierr /= 0 )
return
113 allocate (this%lat(this%nlat))
114 allocate (this%pres(this%nlev))
115 allocate (this%po3(this%nlev))
116 allocate (this%time(this%ntime+1))
117 allocate (this%data(this%nlat,this%nlev,this%ncf,this%ntime))
119 allocate(lat4(this%nlat), pres4(this%nlev), time4(this%ntime+1))
120 read (fileid, iostat=ierr, iomsg=err_message) this%ncf, this%nlat, this%nlev, this%ntime, lat4, pres4, time4
121 if (ierr /= 0 )
return
124 this%pres(:) = pres4(:)
125 this%po3(:) = log(100.0*this%pres(:))
126 this%lat(:) = lat4(:)
127 this%time(:) = time4(:)
128 deallocate(lat4, pres4, time4)
130 allocate(tempin(this%nlat))
134 read(fileid, iostat=ierr, iomsg=err_message) tempin
135 if (ierr /= 0 )
return
136 this%data(:,i3,i2,i1) = tempin(:)
174 subroutine update_o3prog(this, idx1, idx2, idxh, rjday, idxt1, idxt2, ozpl)
176 integer,
intent(in) :: idx1(:), idx2(:)
177 real(kind_phys),
intent(in) :: idxh(:)
178 real(kind_phys),
intent(in) :: rjday
179 integer,
intent(in) :: idxt1, idxt2
180 real(kind_phys),
intent(out) :: ozpl(:,:,:)
181 integer :: nc, l, j, j1, j2
182 real(kind_phys) :: tem, tx1, tx2
184 tx1 = (this%time(idxt2) - rjday) / (this%time(idxt2) - this%time(idxt1))
189 do j=1,
size(ozpl(:,1,1))
193 ozpl(j,l,nc) = tx1*(tem*this%data(j1,l,nc,idxt1)+idxh(j)*this%data(j2,l,nc,idxt1)) &
194 + tx2*(tem*this%data(j1,l,nc,idxt2)+idxh(j)*this%data(j2,l,nc,idxt2))
202 subroutine run_o3prog_2015(this, con_1ovg, dt, p, t, dp, ozpl, oz, do_diag, do3_dt_prd, &
203 do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz)
205 real(kind_phys),
intent(in) :: &
207 real(kind_phys),
intent(in) :: &
209 real(kind_phys),
intent(in),
dimension(:,:) :: &
213 real(kind_phys),
intent(in),
dimension(:,:,:) :: &
215 real(kind_phys),
intent(inout),
dimension(:,:) :: &
217 logical,
intent(in) :: do_diag
218 real(kind_phys),
intent(inout),
dimension(:,:),
optional :: &
224 integer :: k, kmax, kmin, iLev, iCol, nCol, nLev, iCf
225 logical,
dimension(size(p,1)) :: flg
226 real(kind_phys) :: pmax, pmin, tem, temp
227 real(kind_phys),
dimension(size(p,1)) :: wk1, wk2, wk3, ozib
228 real(kind_phys),
dimension(size(p,1),this%ncf) :: prod
229 real(kind_phys),
dimension(size(p,1),size(p,2)) :: ozi
230 real(kind_phys),
dimension(size(p,1),size(p,2)+1) :: colo3, coloz
239 colo3(:,nlev+1) = 0.0
240 coloz(:,nlev+1) = 0.0
247 wk1(icol) = log(p(icol,ilev))
248 pmin = min(wk1(icol), pmin)
249 pmax = max(wk1(icol), pmax)
250 prod(icol,:) = 0._kind_phys
255 if (pmin < this%po3(k)) kmax = k
256 if (pmax < this%po3(k)) kmin = k
260 temp = 1.0 / (this%po3(k) - this%po3(k+1))
263 if (wk1(icol) < this%po3(k) .and. wk1(icol) >= this%po3(k+1))
then
265 wk2(icol) = (wk1(icol) - this%po3(k+1)) * temp
266 wk3(icol) = 1.0 - wk2(icol)
272 prod(icol,icf) = wk2(icol) * ozpl(icol,k,icf) + wk3(icol) * ozpl(icol,k+1,icf)
280 if (wk1(icol) < this%po3(this%nlev))
then
281 prod(icol,icf) = ozpl(icol,this%nlev,icf)
283 if (wk1(icol) >= this%po3(1))
then
284 prod(icol,icf) = ozpl(icol,1,icf)
289 colo3(icol,ilev) = colo3(icol,ilev+1) + ozi(icol,ilev) * dp(icol,ilev)*con_1ovg
290 coloz(icol,ilev) = coloz(icol,ilev+1) + prod(icol,6) * dp(icol,ilev)*con_1ovg
291 prod(icol,2) = min(prod(icol,2), 0.0)
294 ozib(icol) = ozi(icol,ilev)
295 tem = prod(icol,1) - prod(icol,2) * prod(icol,6) &
296 + prod(icol,3) * (t(icol,ilev) - prod(icol,5)) &
297 + prod(icol,4) * (colo3(icol,ilev)-coloz(icol,ilev))
298 oz(icol,ilev) = (ozib(icol) + tem*dt) / (1.0 - prod(icol,2)*dt)
303 do3_dt_prd(:,ilev) = prod(:,1) * dt
304 do3_dt_ozmx(:,ilev) = prod(:,2) * (oz(:,ilev) - prod(:,6)) * dt
305 do3_dt_temp(:,ilev) = prod(:,3)*(t(:,ilev)-prod(:,5))*dt
306 do3_dt_ohoz(:,ilev) = prod(:,4) * (colo3(:,ilev)-coloz(:,ilev))*dt
314 subroutine run_o3prog_2006(this, con_1ovg, dt, p, t, dp, ozpl, oz, do_diag, do3_dt_prd, &
315 do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz)
317 real(kind_phys),
intent(in) :: &
319 real(kind_phys),
intent(in) :: &
321 real(kind_phys),
intent(in),
dimension(:,:) :: &
325 real(kind_phys),
intent(in),
dimension(:,:,:) :: &
327 real(kind_phys),
intent(inout),
dimension(:,:) :: &
329 logical,
intent(in) :: do_diag
330 real(kind_phys),
intent(inout),
dimension(:,:) :: &
337 integer :: k, kmax, kmin, iLev, iCol, nCol, nLev, iCf
338 logical,
dimension(size(p,1)) :: flg
339 real(kind_phys) :: pmax, pmin, tem, temp
340 real(kind_phys),
dimension(size(p,1)) :: wk1, wk2, wk3, ozib
341 real(kind_phys),
dimension(size(p,1),this%ncf) :: prod
342 real(kind_phys),
dimension(size(p,1),size(p,2)) :: ozi
343 real(kind_phys),
dimension(size(p,1),size(p,2)+1) :: colo3, coloz
353 if (this%ncf > 2)
then
354 colo3(:,nlev+1) = 0.0
357 colo3(icol,ilev) = colo3(icol,ilev+1) + ozi(icol,ilev) * dp(icol,ilev) * con_1ovg
368 wk1(icol) = log(p(icol,ilev))
369 pmin = min(wk1(icol), pmin)
370 pmax = max(wk1(icol), pmax)
371 prod(icol,:) = 0._kind_phys
376 if (pmin < this%po3(k)) kmax = k
377 if (pmax < this%po3(k)) kmin = k
381 temp = 1.0 / (this%po3(k) - this%po3(k+1))
384 if (wk1(icol) < this%po3(k) .and. wk1(icol) >= this%po3(k+1))
then
386 wk2(icol) = (wk1(icol) - this%po3(k+1)) * temp
387 wk3(icol) = 1.0 - wk2(icol)
393 prod(icol,icf) = wk2(icol) * ozpl(icol,k,icf) + wk3(icol) * ozpl(icol,k+1,icf)
401 if (wk1(icol) < this%po3(this%nlev))
then
402 prod(icol,icf) = ozpl(icol,this%nlev,icf)
404 if (wk1(icol) >= this%po3(1))
then
405 prod(icol,icf) = ozpl(icol,1,icf)
410 if (this%ncf == 2)
then
412 ozib(icol) = ozi(icol,ilev)
413 oz(icol,ilev) = (ozib(icol) + prod(icol,1)*dt) / (1.0 + prod(icol,2)*dt)
417 if (this%ncf == 4)
then
419 ozib(icol) = ozi(icol,ilev)
420 tem = prod(icol,1) + prod(icol,3)*t(icol,ilev) + prod(icol,4)*colo3(icol,ilev+1)
421 oz(icol,ilev) = (ozib(icol) + tem*dt) / (1.0 + prod(icol,2)*dt)
427 do3_dt_prd(:,ilev) = prod(:,1)*dt
428 do3_dt_ozmx(:,ilev) = (oz(:,ilev) - ozib(:))
429 do3_dt_temp(:,ilev) = prod(:,3) * t(:,ilev) * dt
430 do3_dt_ohoz(:,ilev) = prod(:,4) * colo3(:,ilev) * dt
438 subroutine run_o3clim(this, lat, prslk, con_pi, oz)
440 real(kind_phys),
intent(in) :: &
442 real(kind_phys),
intent(in),
dimension(:) :: &
444 real(kind_phys),
intent(in),
dimension(:,:) :: &
446 real(kind_phys),
intent(out),
dimension(:,:) :: &
449 integer :: nCol, iCol, nLev, iLev, j, j1, j2, l, ll
450 real(kind_phys) :: elte, deglat, tem, tem1, tem2, tem3, tem4, temp
451 real(kind_phys),
allocatable :: o3i(:,:),wk1(:)
454 ncol =
size(prslk(:,1))
455 nlev =
size(prslk(1,:))
456 allocate(o3i(ncol, this%nlevc),wk1(ncol))
459 top_at_1 = (prslk(1,1) .lt. prslk(1, nlev))
461 elte = this%blatc + (this%nlatc-1)*this%dphiozc
464 deglat = lat(icol) * 180.0 / con_pi
465 if (deglat > this%blatc .and. deglat < elte)
then
466 tem1 = (deglat - this%blatc) / this%dphiozc + 1
470 elseif (deglat <= this%blatc)
then
474 elseif (deglat >= elte)
then
482 tem3 = tem2*this%datac(j1,j,this%k1oz) + tem1*this%datac(j2,j,this%k1oz)
483 tem4 = tem2*this%datac(j1,j,this%k2oz) + tem1*this%datac(j2,j,this%k2oz)
484 o3i(icol,j) = tem4*this%facoz + tem3*(1.0 - this%facoz)
490 if (.not. top_at_1) ll = nlev - ilev + 1
493 wk1(icol) = prslk(icol,ll)
496 do j = 1, this%nlevc-1
497 temp = 1.0 / (this%pkstr(j+1) - this%pkstr(j))
499 if (wk1(icol) > this%pkstr(j) .and. wk1(icol) <= this%pkstr(j+1))
then
500 tem = (this%pkstr(j+1) - wk1(icol)) * temp
501 oz(icol,ll) = tem * o3i(icol,j) + (1.0 - tem) * o3i(icol,j+1)
507 if (wk1(icol) > this%pkstr(this%nlevc)) oz(icol,ll) = o3i(icol,this%nlevc)
508 if (wk1(icol) < this%pkstr(1)) oz(icol,ll) = o3i(icol,1)
516 function load_o3clim(this, file, fileID)
result (err_message)
518 integer,
intent(in) :: fileid
519 character(len=*),
intent(in) :: file
520 character(len=128) :: err_message
523 real(kind=4) :: blatc4
524 integer :: ilev, ilat, imo, ierr
525 real(kind=4), allocatable :: o3clim4(:,:,:), pstr4(:)
526 integer,
allocatable :: imond(:), ilatt(:,:)
531 open(unit=fileid,file=trim(file),form=
'unformatted',convert=
'big_endian', iostat=ierr, iomsg=err_message)
532 if (ierr /= 0 )
return
533 read (fileid,
end=101,iostat=ierr,iomsg=err_message) this%nlatc, this%nlevc, this%ntimec, blatc4
534 if (ierr /= 0 )
return
536101
if (this%nlevc < 10 .or. this%nlevc > 100)
then
548 this%dphiozc = -(this%blatc+this%blatc)/(this%nlatc-1)
550 allocate (o3clim4(this%nlatc,this%nlevc,12), pstr4(this%nlevc), imond(12), ilatt(this%nlatc,12) )
552 allocate (this%pkstr(this%nlevc), this%pstr(this%nlevc), this%datac(this%nlatc,this%nlevc,12))
553 if ( this%nlevc == 17 )
then
554 do ilev = 1, this%nlevc
555 read (fileid,15,iostat=ierr,iomsg=err_message) pstr4(ilev)
556 if (ierr /= 0 )
return
561 do ilat = 1, this%nlatc
562 read (fileid,16,iostat=ierr,iomsg=err_message) imond(imo), ilatt(ilat,imo), (o3clim4(ilat,ilev,imo),ilev=1,10)
563 if (ierr /= 0 )
return
56416
format(i2,i4,10f6.2)
565 read (fileid,20,iostat=ierr,iomsg=err_message) (o3clim4(ilat,ilev,imo),ilev=11,this%nlevc)
566 if (ierr /= 0 )
return
572 do ilev = 1, this%nlevc
573 read (fileid) pstr4(ilev)
577 do ilev = 1, this%nlevc
578 read (fileid,iostat=ierr,iomsg=err_message) (o3clim4(ilat,ilev,imo),ilat=1,this%nlatc)
579 if (ierr /= 0 )
return
585 do ilev = 1, this%nlevc
586 do ilat = 1, this%nlatc
587 this%datac(ilat,ilev,imo) = o3clim4(ilat,ilev,imo) * 1.655e-6
592 do ilev = 1, this%nlevc
593 this%pstr(ilev) = pstr4(ilev)
594 this%pkstr(ilev) = fpkapx(this%pstr(ilev)*100.0)
601 subroutine update_o3clim(this, imon, iday, ihour, loz1st)
603 integer,
intent(in) :: imon, iday, ihour
604 logical,
intent(in) :: loz1st
606 integer :: midmon=15, midm=15, midp=45, id
607 integer,
parameter,
dimension(13) :: mdays = (/31,28,31,30,31,30,31,31,30,31,30,31,30/)
610 midmon = mdays(imon)/2 + 1
611 change = loz1st .or. ( (iday==midmon) .and. (ihour==0) )
614 if ( iday < midmon )
then
615 this%k1oz = mod(imon+10, 12) + 1
616 midm = mdays(this%k1oz)/2 + 1
618 midp = mdays(this%k1oz) + midmon
622 this%k2oz = mod(imon, 12) + 1
623 midp = mdays(this%k2oz)/2 + 1 + mdays(this%k1oz)
627 if (iday < midmon)
then
628 id = iday + mdays(this%k1oz)
633 this%facoz = float(id - midm) / float(midp - midm)