88 function load_o3prog(this, file, fileID)
result (err_message)
90 integer,
intent(in) :: fileid
91 character(len=*),
intent(in) :: file
92 character(len=128) :: err_message
93 integer :: i1, i2, i3, ierr
94 real(kind=4), dimension(:),
allocatable :: lat4, pres4, time4, tempin
95 real(kind=4) :: blatc4
101 open(unit=fileid,file=trim(file), form=
'unformatted', convert=
'big_endian', iostat=ierr, iomsg=err_message)
102 if (ierr /= 0 )
return
103 read (fileid, iostat=ierr, iomsg=err_message) this%ncf, this%nlat, this%nlev, this%ntime
104 if (ierr /= 0 )
return
107 allocate (this%lat(this%nlat))
108 allocate (this%pres(this%nlev))
109 allocate (this%po3(this%nlev))
110 allocate (this%time(this%ntime+1))
111 allocate (this%data(this%nlat,this%nlev,this%ncf,this%ntime))
113 allocate(lat4(this%nlat), pres4(this%nlev), time4(this%ntime+1))
114 read (fileid, iostat=ierr, iomsg=err_message) this%ncf, this%nlat, this%nlev, this%ntime, lat4, pres4, time4
115 if (ierr /= 0 )
return
118 this%pres(:) = pres4(:)
119 this%po3(:) = log(100.0*this%pres(:))
120 this%lat(:) = lat4(:)
121 this%time(:) = time4(:)
122 deallocate(lat4, pres4, time4)
124 allocate(tempin(this%nlat))
128 read(fileid, iostat=ierr, iomsg=err_message) tempin
129 if (ierr /= 0 )
return
130 this%data(:,i3,i2,i1) = tempin(:)
168 subroutine update_o3prog(this, idx1, idx2, idxh, rjday, idxt1, idxt2, ozpl)
170 integer,
intent(in) :: idx1(:), idx2(:)
171 real(kind_phys),
intent(in) :: idxh(:)
172 real(kind_phys),
intent(in) :: rjday
173 integer,
intent(in) :: idxt1, idxt2
174 real(kind_phys),
intent(out) :: ozpl(:,:,:)
175 integer :: nc, l, j, j1, j2
176 real(kind_phys) :: tem, tx1, tx2
178 tx1 = (this%time(idxt2) - rjday) / (this%time(idxt2) - this%time(idxt1))
183 do j=1,
size(ozpl(:,1,1))
187 ozpl(j,l,nc) = tx1*(tem*this%data(j1,l,nc,idxt1)+idxh(j)*this%data(j2,l,nc,idxt1)) &
188 + tx2*(tem*this%data(j1,l,nc,idxt2)+idxh(j)*this%data(j2,l,nc,idxt2))
196 subroutine run_o3prog_2015(this, con_1ovg, dt, p, t, dp, ozpl, oz, do_diag, do3_dt_prd, &
197 do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz)
199 real(kind_phys),
intent(in) :: &
201 real(kind_phys),
intent(in) :: &
203 real(kind_phys),
intent(in),
dimension(:,:) :: &
207 real(kind_phys),
intent(in),
dimension(:,:,:) :: &
209 real(kind_phys),
intent(inout),
dimension(:,:) :: &
211 logical,
intent(in) :: do_diag
212 real(kind_phys),
intent(inout),
dimension(:,:),
optional :: &
218 integer :: k, kmax, kmin, iLev, iCol, nCol, nLev, iCf
219 logical,
dimension(size(p,1)) :: flg
220 real(kind_phys) :: pmax, pmin, tem, temp
221 real(kind_phys),
dimension(size(p,1)) :: wk1, wk2, wk3, ozib
222 real(kind_phys),
dimension(size(p,1),this%ncf) :: prod
223 real(kind_phys),
dimension(size(p,1),size(p,2)) :: ozi
224 real(kind_phys),
dimension(size(p,1),size(p,2)+1) :: colo3, coloz
233 colo3(:,nlev+1) = 0.0
234 coloz(:,nlev+1) = 0.0
241 wk1(icol) = log(p(icol,ilev))
242 pmin = min(wk1(icol), pmin)
243 pmax = max(wk1(icol), pmax)
244 prod(icol,:) = 0._kind_phys
249 if (pmin < this%po3(k)) kmax = k
250 if (pmax < this%po3(k)) kmin = k
254 temp = 1.0 / (this%po3(k) - this%po3(k+1))
257 if (wk1(icol) < this%po3(k) .and. wk1(icol) >= this%po3(k+1))
then
259 wk2(icol) = (wk1(icol) - this%po3(k+1)) * temp
260 wk3(icol) = 1.0 - wk2(icol)
266 prod(icol,icf) = wk2(icol) * ozpl(icol,k,icf) + wk3(icol) * ozpl(icol,k+1,icf)
274 if (wk1(icol) < this%po3(this%nlev))
then
275 prod(icol,icf) = ozpl(icol,this%nlev,icf)
277 if (wk1(icol) >= this%po3(1))
then
278 prod(icol,icf) = ozpl(icol,1,icf)
283 colo3(icol,ilev) = colo3(icol,ilev+1) + ozi(icol,ilev) * dp(icol,ilev)*con_1ovg
284 coloz(icol,ilev) = coloz(icol,ilev+1) + prod(icol,6) * dp(icol,ilev)*con_1ovg
285 prod(icol,2) = min(prod(icol,2), 0.0)
288 ozib(icol) = ozi(icol,ilev)
289 tem = prod(icol,1) - prod(icol,2) * prod(icol,6) &
290 + prod(icol,3) * (t(icol,ilev) - prod(icol,5)) &
291 + prod(icol,4) * (colo3(icol,ilev)-coloz(icol,ilev))
292 oz(icol,ilev) = (ozib(icol) + tem*dt) / (1.0 - prod(icol,2)*dt)
297 do3_dt_prd(:,ilev) = prod(:,1) * dt
298 do3_dt_ozmx(:,ilev) = prod(:,2) * (oz(:,ilev) - prod(:,6)) * dt
299 do3_dt_temp(:,ilev) = prod(:,3)*(t(:,ilev)-prod(:,5))*dt
300 do3_dt_ohoz(:,ilev) = prod(:,4) * (colo3(:,ilev)-coloz(:,ilev))*dt
308 subroutine run_o3prog_2006(this, con_1ovg, dt, p, t, dp, ozpl, oz, do_diag, do3_dt_prd, &
309 do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz)
311 real(kind_phys),
intent(in) :: &
313 real(kind_phys),
intent(in) :: &
315 real(kind_phys),
intent(in),
dimension(:,:) :: &
319 real(kind_phys),
intent(in),
dimension(:,:,:) :: &
321 real(kind_phys),
intent(inout),
dimension(:,:) :: &
323 logical,
intent(in) :: do_diag
324 real(kind_phys),
intent(inout),
dimension(:,:) :: &
331 integer :: k, kmax, kmin, iLev, iCol, nCol, nLev, iCf
332 logical,
dimension(size(p,1)) :: flg
333 real(kind_phys) :: pmax, pmin, tem, temp
334 real(kind_phys),
dimension(size(p,1)) :: wk1, wk2, wk3, ozib
335 real(kind_phys),
dimension(size(p,1),this%ncf) :: prod
336 real(kind_phys),
dimension(size(p,1),size(p,2)) :: ozi
337 real(kind_phys),
dimension(size(p,1),size(p,2)+1) :: colo3, coloz
347 if (this%ncf > 2)
then
348 colo3(:,nlev+1) = 0.0
351 colo3(icol,ilev) = colo3(icol,ilev+1) + ozi(icol,ilev) * dp(icol,ilev) * con_1ovg
362 wk1(icol) = log(p(icol,ilev))
363 pmin = min(wk1(icol), pmin)
364 pmax = max(wk1(icol), pmax)
365 prod(icol,:) = 0._kind_phys
370 if (pmin < this%po3(k)) kmax = k
371 if (pmax < this%po3(k)) kmin = k
375 temp = 1.0 / (this%po3(k) - this%po3(k+1))
378 if (wk1(icol) < this%po3(k) .and. wk1(icol) >= this%po3(k+1))
then
380 wk2(icol) = (wk1(icol) - this%po3(k+1)) * temp
381 wk3(icol) = 1.0 - wk2(icol)
387 prod(icol,icf) = wk2(icol) * ozpl(icol,k,icf) + wk3(icol) * ozpl(icol,k+1,icf)
395 if (wk1(icol) < this%po3(this%nlev))
then
396 prod(icol,icf) = ozpl(icol,this%nlev,icf)
398 if (wk1(icol) >= this%po3(1))
then
399 prod(icol,icf) = ozpl(icol,1,icf)
404 if (this%ncf == 2)
then
406 ozib(icol) = ozi(icol,ilev)
407 oz(icol,ilev) = (ozib(icol) + prod(icol,1)*dt) / (1.0 + prod(icol,2)*dt)
411 if (this%ncf == 4)
then
413 ozib(icol) = ozi(icol,ilev)
414 tem = prod(icol,1) + prod(icol,3)*t(icol,ilev) + prod(icol,4)*colo3(icol,ilev+1)
415 oz(icol,ilev) = (ozib(icol) + tem*dt) / (1.0 + prod(icol,2)*dt)
421 do3_dt_prd(:,ilev) = prod(:,1)*dt
422 do3_dt_ozmx(:,ilev) = (oz(:,ilev) - ozib(:))
423 do3_dt_temp(:,ilev) = prod(:,3) * t(:,ilev) * dt
424 do3_dt_ohoz(:,ilev) = prod(:,4) * colo3(:,ilev) * dt
432 subroutine run_o3clim(this, lat, prslk, con_pi, oz)
434 real(kind_phys),
intent(in) :: &
436 real(kind_phys),
intent(in),
dimension(:) :: &
438 real(kind_phys),
intent(in),
dimension(:,:) :: &
440 real(kind_phys),
intent(out),
dimension(:,:) :: &
443 integer :: nCol, iCol, nLev, iLev, j, j1, j2, l, ll
444 real(kind_phys) :: elte, deglat, tem, tem1, tem2, tem3, tem4, temp
445 real(kind_phys),
allocatable :: o3i(:,:),wk1(:)
448 ncol =
size(prslk(:,1))
449 nlev =
size(prslk(1,:))
450 allocate(o3i(ncol, this%nlevc),wk1(ncol))
453 top_at_1 = (prslk(1,1) .lt. prslk(1, nlev))
455 elte = this%blatc + (this%nlatc-1)*this%dphiozc
458 deglat = lat(icol) * 180.0 / con_pi
459 if (deglat > this%blatc .and. deglat < elte)
then
460 tem1 = (deglat - this%blatc) / this%dphiozc + 1
464 elseif (deglat <= this%blatc)
then
468 elseif (deglat >= elte)
then
476 tem3 = tem2*this%datac(j1,j,this%k1oz) + tem1*this%datac(j2,j,this%k1oz)
477 tem4 = tem2*this%datac(j1,j,this%k2oz) + tem1*this%datac(j2,j,this%k2oz)
478 o3i(icol,j) = tem4*this%facoz + tem3*(1.0 - this%facoz)
484 if (.not. top_at_1) ll = nlev - ilev + 1
487 wk1(icol) = prslk(icol,ll)
490 do j = 1, this%nlevc-1
491 temp = 1.0 / (this%pkstr(j+1) - this%pkstr(j))
493 if (wk1(icol) > this%pkstr(j) .and. wk1(icol) <= this%pkstr(j+1))
then
494 tem = (this%pkstr(j+1) - wk1(icol)) * temp
495 oz(icol,ll) = tem * o3i(icol,j) + (1.0 - tem) * o3i(icol,j+1)
501 if (wk1(icol) > this%pkstr(this%nlevc)) oz(icol,ll) = o3i(icol,this%nlevc)
502 if (wk1(icol) < this%pkstr(1)) oz(icol,ll) = o3i(icol,1)
510 function load_o3clim(this, file, fileID)
result (err_message)
512 integer,
intent(in) :: fileid
513 character(len=*),
intent(in) :: file
514 character(len=128) :: err_message
517 real(kind=4) :: blatc4
518 integer :: ilev, ilat, imo, ierr
519 real(kind=4), allocatable :: o3clim4(:,:,:), pstr4(:)
520 integer,
allocatable :: imond(:), ilatt(:,:)
525 open(unit=fileid,file=trim(file),form=
'unformatted',convert=
'big_endian', iostat=ierr, iomsg=err_message)
526 if (ierr /= 0 )
return
527 read (fileid,
end=101,iostat=ierr,iomsg=err_message) this%nlatc, this%nlevc, this%ntimec, blatc4
528 if (ierr /= 0 )
return
530101
if (this%nlevc < 10 .or. this%nlevc > 100)
then
542 this%dphiozc = -(this%blatc+this%blatc)/(this%nlatc-1)
544 allocate (o3clim4(this%nlatc,this%nlevc,12), pstr4(this%nlevc), imond(12), ilatt(this%nlatc,12) )
546 allocate (this%pkstr(this%nlevc), this%pstr(this%nlevc), this%datac(this%nlatc,this%nlevc,12))
547 if ( this%nlevc == 17 )
then
548 do ilev = 1, this%nlevc
549 read (fileid,15,iostat=ierr,iomsg=err_message) pstr4(ilev)
550 if (ierr /= 0 )
return
555 do ilat = 1, this%nlatc
556 read (fileid,16,iostat=ierr,iomsg=err_message) imond(imo), ilatt(ilat,imo), (o3clim4(ilat,ilev,imo),ilev=1,10)
557 if (ierr /= 0 )
return
55816
format(i2,i4,10f6.2)
559 read (fileid,20,iostat=ierr,iomsg=err_message) (o3clim4(ilat,ilev,imo),ilev=11,this%nlevc)
560 if (ierr /= 0 )
return
566 do ilev = 1, this%nlevc
567 read (fileid) pstr4(ilev)
571 do ilev = 1, this%nlevc
572 read (fileid,iostat=ierr,iomsg=err_message) (o3clim4(ilat,ilev,imo),ilat=1,this%nlatc)
573 if (ierr /= 0 )
return
579 do ilev = 1, this%nlevc
580 do ilat = 1, this%nlatc
581 this%datac(ilat,ilev,imo) = o3clim4(ilat,ilev,imo) * 1.655e-6
586 do ilev = 1, this%nlevc
587 this%pstr(ilev) = pstr4(ilev)
588 this%pkstr(ilev) = fpkapx(this%pstr(ilev)*100.0)
595 subroutine update_o3clim(this, imon, iday, ihour, loz1st)
597 integer,
intent(in) :: imon, iday, ihour
598 logical,
intent(in) :: loz1st
600 integer :: midmon=15, midm=15, midp=45, id
601 integer,
parameter,
dimension(13) :: mdays = (/31,28,31,30,31,30,31,31,30,31,30,31,30/)
604 midmon = mdays(imon)/2 + 1
605 change = loz1st .or. ( (iday==midmon) .and. (ihour==0) )
608 if ( iday < midmon )
then
609 this%k1oz = mod(imon+10, 12) + 1
610 midm = mdays(this%k1oz)/2 + 1
612 midp = mdays(this%k1oz) + midmon
616 this%k2oz = mod(imon, 12) + 1
617 midp = mdays(this%k2oz)/2 + 1 + mdays(this%k1oz)
621 if (iday < midmon)
then
622 id = iday + mdays(this%k1oz)
627 this%facoz = float(id - midm) / float(midp - midm)