CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_ozphys.F90
1
3! #########################################################################################
7!
40! #########################################################################################
42 use machine, only : kind_phys
43 use funcphys, only : fpkapx
44 implicit none
45
46 public ty_ozphys
47
48! #########################################################################################
55! #########################################################################################
57 ! Prognostic ozone.
58 integer :: nlat
59 integer :: nlev
60 integer :: ntime
61 integer :: ncf
62 real(kind_phys), allocatable :: lat(:)
63 real(kind_phys), allocatable :: pres(:)
64 real(kind_phys), allocatable :: po3(:)
65 real(kind_phys), allocatable :: time(:)
66 real(kind_phys), allocatable :: data(:,:,:,:)
67 ! Climotological ozone.
68 integer :: nlatc
69 integer :: nlevc
70 integer :: ntimec
71 real(kind_phys) :: blatc
72 real(kind_phys) :: dphiozc
73 real(kind_phys), allocatable :: pkstr(:)
74 real(kind_phys), allocatable :: pstr(:)
75 real(kind_phys), allocatable :: datac(:,:,:)
76 integer :: k1oz
77 integer :: k2oz
78 real(kind_phys) :: facoz
79 contains
80 procedure, public :: load_o3prog
81 procedure, public :: setup_o3prog
82 procedure, public :: update_o3prog
83 procedure, public :: run_o3prog_2015
84 procedure, public :: run_o3prog_2006
85 !
86 procedure, public :: load_o3clim
87 procedure, public :: update_o3clim
88 procedure, public :: run_o3clim
89 end type ty_ozphys
90
91contains
92
94 function load_o3prog(this, file, fileID) result (err_message)
95 class(ty_ozphys), intent(inout) :: this
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
102
103 ! initialize error message
104 err_message = ""
105
106 ! Get dimensions from data file
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
111 rewind(fileid)
112
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))
118
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
122
123 ! Store
124 this%pres(:) = pres4(:)
125 this%po3(:) = log(100.0*this%pres(:)) ! from mb to ln(Pa)
126 this%lat(:) = lat4(:)
127 this%time(:) = time4(:)
128 deallocate(lat4, pres4, time4)
129
130 allocate(tempin(this%nlat))
131 do i1=1,this%ntime
132 do i2=1,this%ncf
133 do i3=1,this%nlev
134 read(fileid, iostat=ierr, iomsg=err_message) tempin
135 if (ierr /= 0 ) return
136 this%data(:,i3,i2,i1) = tempin(:)
137 enddo
138 enddo
139 enddo
140 deallocate(tempin)
141 close(fileid)
142
143 end function load_o3prog
144
147 subroutine setup_o3prog(this, lat, idx1, idx2, idxh)
148 class(ty_ozphys), intent(in) :: this
149 real(kind_phys), intent(in) :: lat(:)
150 integer, intent(out) :: idx1(:), idx2(:)
151 real(kind_phys), intent(out) :: idxh(:)
152 integer :: i,j
153
154 do j=1,size(lat)
155 idx2(j) = this%nlat + 1
156 do i=1,this%nlat
157 if (lat(j) < this%lat(i)) then
158 idx2(j) = i
159 exit
160 endif
161 enddo
162 idx1(j) = max(idx2(j)-1,1)
163 idx2(j) = min(idx2(j),this%nlat)
164 if (idx2(j) .ne. idx1(j)) then
165 idxh(j) = (lat(j) - this%lat(idx1(j))) / (this%lat(idx2(j)) - this%lat(idx1(j)))
166 else
167 idxh(j) = 1.0
168 endif
169 enddo
170
171 end subroutine setup_o3prog
172
174 subroutine update_o3prog(this, idx1, idx2, idxh, rjday, idxt1, idxt2, ozpl)
175 class(ty_ozphys), intent(in) :: this
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
183
184 tx1 = (this%time(idxt2) - rjday) / (this%time(idxt2) - this%time(idxt1))
185 tx2 = 1.0 - tx1
186
187 do nc=1,this%ncf
188 do l=1,this%nlev
189 do j=1,size(ozpl(:,1,1))
190 j1 = idx1(j)
191 j2 = idx2(j)
192 tem = 1.0 - idxh(j)
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))
195 enddo
196 enddo
197 enddo
198
199 end subroutine update_o3prog
200
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)
204 class(ty_ozphys), intent(in) :: this
205 real(kind_phys), intent(in) :: &
206 con_1ovg ! Physical constant: One divided by gravitational acceleration (m-1 s2)
207 real(kind_phys), intent(in) :: &
208 dt ! Model timestep (sec)
209 real(kind_phys), intent(in), dimension(:,:) :: &
210 p, & ! Model Pressure (Pa)
211 t, & ! Model temperature (K)
212 dp ! Model layer thickness (Pa)
213 real(kind_phys), intent(in), dimension(:,:,:) :: &
214 ozpl ! Ozone forcing data
215 real(kind_phys), intent(inout), dimension(:,:) :: &
216 oz ! Ozone concentration updated by physics
217 logical, intent(in) :: do_diag
218 real(kind_phys), intent(inout), dimension(:,:), optional :: &
219 do3_dt_prd, & ! Physics tendency: production and loss effect
220 do3_dt_ozmx, & ! Physics tendency: ozone mixing ratio effect
221 do3_dt_temp, & ! Physics tendency: temperature effect
222 do3_dt_ohoz ! Physics tendency: overhead ozone effect
223
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
231
232 ! Dimensions
233 ncol = size(p,1)
234 nlev = size(p,2)
235
236 ! Temporaries
237 ozi = oz
238
239 colo3(:,nlev+1) = 0.0
240 coloz(:,nlev+1) = 0.0
241
242 do ilev=nlev,1,-1
243 pmin = 1.0e10
244 pmax = -1.0e10
245
246 do icol=1,ncol
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
251 enddo
252 kmax = 1
253 kmin = 1
254 do k=1,this%nlev-1
255 if (pmin < this%po3(k)) kmax = k
256 if (pmax < this%po3(k)) kmin = k
257 enddo
258 !
259 do k=kmin,kmax
260 temp = 1.0 / (this%po3(k) - this%po3(k+1))
261 do icol=1,ncol
262 flg(icol) = .false.
263 if (wk1(icol) < this%po3(k) .and. wk1(icol) >= this%po3(k+1)) then
264 flg(icol) = .true.
265 wk2(icol) = (wk1(icol) - this%po3(k+1)) * temp
266 wk3(icol) = 1.0 - wk2(icol)
267 endif
268 enddo
269 do icf=1,this%ncf
270 do icol=1,ncol
271 if (flg(icol)) then
272 prod(icol,icf) = wk2(icol) * ozpl(icol,k,icf) + wk3(icol) * ozpl(icol,k+1,icf)
273 endif
274 enddo
275 enddo
276 enddo
277
278 do icf=1,this%ncf
279 do icol=1,ncol
280 if (wk1(icol) < this%po3(this%nlev)) then
281 prod(icol,icf) = ozpl(icol,this%nlev,icf)
282 endif
283 if (wk1(icol) >= this%po3(1)) then
284 prod(icol,icf) = ozpl(icol,1,icf)
285 endif
286 enddo
287 enddo
288 do icol=1,ncol
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)
292 enddo
293 do icol=1,ncol
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)
299 enddo
300
301 ! Diagnostics (optional)
302 if (do_diag) then
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
307 endif
308 enddo
309
310 return
311 end subroutine run_o3prog_2015
312
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)
316 class(ty_ozphys), intent(in) :: this
317 real(kind_phys), intent(in) :: &
318 con_1ovg ! Physical constant: One divided by gravitational acceleration (m-1 s2)
319 real(kind_phys), intent(in) :: &
320 dt ! Model timestep (sec)
321 real(kind_phys), intent(in), dimension(:,:) :: &
322 p, & ! Model Pressure (Pa)
323 t, & ! Model temperature (K)
324 dp ! Model layer thickness (Pa)
325 real(kind_phys), intent(in), dimension(:,:,:) :: &
326 ozpl ! Ozone forcing data
327 real(kind_phys), intent(inout), dimension(:,:) :: &
328 oz ! Ozone concentration updated by physics
329 logical, intent(in) :: do_diag
330 real(kind_phys), intent(inout), dimension(:,:) :: &
331 do3_dt_prd, & ! Physics tendency: production and loss effect
332 do3_dt_ozmx, & ! Physics tendency: ozone mixing ratio effect
333 do3_dt_temp, & ! Physics tendency: temperature effect
334 do3_dt_ohoz ! Physics tendency: overhead ozone effect
335
336 ! Locals
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
344
345 ! Dimensions
346 ncol = size(p,1)
347 nlev = size(p,2)
348
349 ! Temporaries
350 ozi = oz
351
353 if (this%ncf > 2) then
354 colo3(:,nlev+1) = 0.0
355 do ilev=nlev,1,-1
356 do icol=1,ncol
357 colo3(icol,ilev) = colo3(icol,ilev+1) + ozi(icol,ilev) * dp(icol,ilev) * con_1ovg
358 enddo
359 enddo
360 endif
361
363 do ilev=1,nlev
364 pmin = 1.0e10
365 pmax = -1.0e10
366
367 do icol=1,ncol
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
372 enddo
373 kmax = 1
374 kmin = 1
375 do k=1,this%nlev-1
376 if (pmin < this%po3(k)) kmax = k
377 if (pmax < this%po3(k)) kmin = k
378 enddo
379
380 do k=kmin,kmax
381 temp = 1.0 / (this%po3(k) - this%po3(k+1))
382 do icol=1,ncol
383 flg(icol) = .false.
384 if (wk1(icol) < this%po3(k) .and. wk1(icol) >= this%po3(k+1)) then
385 flg(icol) = .true.
386 wk2(icol) = (wk1(icol) - this%po3(k+1)) * temp
387 wk3(icol) = 1.0 - wk2(icol)
388 endif
389 enddo
390 do icf=1,this%ncf
391 do icol=1,ncol
392 if (flg(icol)) then
393 prod(icol,icf) = wk2(icol) * ozpl(icol,k,icf) + wk3(icol) * ozpl(icol,k+1,icf)
394 endif
395 enddo
396 enddo
397 enddo
398
399 do icf=1,this%ncf
400 do icol=1,ncol
401 if (wk1(icol) < this%po3(this%nlev)) then
402 prod(icol,icf) = ozpl(icol,this%nlev,icf)
403 endif
404 if (wk1(icol) >= this%po3(1)) then
405 prod(icol,icf) = ozpl(icol,1,icf)
406 endif
407 enddo
408 enddo
409
410 if (this%ncf == 2) then
411 do icol=1,ncol
412 ozib(icol) = ozi(icol,ilev)
413 oz(icol,ilev) = (ozib(icol) + prod(icol,1)*dt) / (1.0 + prod(icol,2)*dt)
414 enddo
415 endif
416
417 if (this%ncf == 4) then
418 do icol=1,ncol
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)
422 enddo
423 endif
424
425 ! Diagnostics (optional)
426 if (do_diag) then
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
431 endif
432 enddo
433
434 return
435 end subroutine run_o3prog_2006
436
438 subroutine run_o3clim(this, lat, prslk, con_pi, oz)
439 class(ty_ozphys), intent(in) :: this
440 real(kind_phys), intent(in) :: &
441 con_pi ! Physics constant: Pi
442 real(kind_phys), intent(in), dimension(:) :: &
443 lat ! Grid latitude
444 real(kind_phys), intent(in), dimension(:,:) :: &
445 prslk ! Exner function
446 real(kind_phys), intent(out), dimension(:,:) :: &
447 oz ! Ozone concentration updated by climotology
448
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(:)
452 logical :: top_at_1
453
454 ncol = size(prslk(:,1))
455 nlev = size(prslk(1,:))
456 allocate(o3i(ncol, this%nlevc),wk1(ncol))
457
458 ! What is vertical ordering?
459 top_at_1 = (prslk(1,1) .lt. prslk(1, nlev))
460
461 elte = this%blatc + (this%nlatc-1)*this%dphiozc
462
463 do icol = 1, ncol
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
467 j1 = tem1
468 j2 = j1 + 1
469 tem1 = tem1 - j1
470 elseif (deglat <= this%blatc) then
471 j1 = 1
472 j2 = 1
473 tem1 = 1.0
474 elseif (deglat >= elte) then
475 j1 = this%nlatc
476 j2 = this%nlatc
477 tem1 = 1.0
478 endif
479
480 tem2 = 1.0 - tem1
481 do j = 1, this%nlevc
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)
485 enddo
486 enddo
487
488 do ilev = 1, nlev
489 ll = ilev
490 if (.not. top_at_1) ll = nlev - ilev + 1
491
492 do icol = 1, ncol
493 wk1(icol) = prslk(icol,ll)
494 enddo
495
496 do j = 1, this%nlevc-1
497 temp = 1.0 / (this%pkstr(j+1) - this%pkstr(j))
498 do icol = 1, ncol
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)
502 endif
503 enddo
504 enddo
505
506 do icol = 1, ncol
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)
509 enddo
510 enddo
511
512 return
513 end subroutine run_o3clim
514
516 function load_o3clim(this, file, fileID) result (err_message)
517 class(ty_ozphys), intent(inout) :: this
518 integer, intent(in) :: fileid
519 character(len=*), intent(in) :: file
520 character(len=128) :: err_message
521
522 ! Locals
523 real(kind=4) :: blatc4
524 integer :: ilev, ilat, imo, ierr
525 real(kind=4), allocatable :: o3clim4(:,:,:), pstr4(:)
526 integer, allocatable :: imond(:), ilatt(:,:)
527
528 ! initialize error message
529 err_message = ""
530
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
535
536101 if (this%nlevc < 10 .or. this%nlevc > 100) then
537 rewind(fileid)
538 this%nlevc = 17
539 this%nlatc = 18
540 this%blatc = -85.0
541 else
542 this%blatc = blatc4
543 endif
544 this%nlat = 2
545 this%nlev = 1
546 this%ntimec = 1
547 this%ncf = 0
548 this%dphiozc = -(this%blatc+this%blatc)/(this%nlatc-1)
549
550 allocate (o3clim4(this%nlatc,this%nlevc,12), pstr4(this%nlevc), imond(12), ilatt(this%nlatc,12) )
551
552 allocate (this%pkstr(this%nlevc), this%pstr(this%nlevc), this%datac(this%nlatc,this%nlevc,12))
553 if ( this%nlevc == 17 ) then ! For the operational ozone climatology
554 do ilev = 1, this%nlevc
555 read (fileid,15,iostat=ierr,iomsg=err_message) pstr4(ilev)
556 if (ierr /= 0 ) return
55715 format(f10.3)
558 enddo
559
560 do imo = 1, 12
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
56720 format(6x,10f6.2)
568 enddo
569 enddo
570 else ! For newer ozone climatology
571 read (fileid)
572 do ilev = 1, this%nlevc
573 read (fileid) pstr4(ilev)
574 enddo
575
576 do imo = 1, 12
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
580 enddo
581 enddo
582 endif ! end if_this%nlevc_block
583
584 do imo = 1, 12
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
588 enddo
589 enddo
590 enddo
591
592 do ilev = 1, this%nlevc
593 this%pstr(ilev) = pstr4(ilev)
594 this%pkstr(ilev) = fpkapx(this%pstr(ilev)*100.0)
595 enddo
596
597 end function load_o3clim
598
601 subroutine update_o3clim(this, imon, iday, ihour, loz1st)
602 class(ty_ozphys), intent(inout) :: this
603 integer, intent(in) :: imon, iday, ihour
604 logical, intent(in) :: loz1st
605
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/)
608 logical :: change
609
610 midmon = mdays(imon)/2 + 1
611 change = loz1st .or. ( (iday==midmon) .and. (ihour==0) )
612
613 if ( change ) then
614 if ( iday < midmon ) then
615 this%k1oz = mod(imon+10, 12) + 1
616 midm = mdays(this%k1oz)/2 + 1
617 this%k2oz = imon
618 midp = mdays(this%k1oz) + midmon
619 else
620 this%k1oz = imon
621 midm = midmon
622 this%k2oz = mod(imon, 12) + 1
623 midp = mdays(this%k2oz)/2 + 1 + mdays(this%k1oz)
624 endif
625 endif
626
627 if (iday < midmon) then
628 id = iday + mdays(this%k1oz)
629 else
630 id = iday
631 endif
632
633 this%facoz = float(id - midm) / float(midp - midm)
634
635 end subroutine update_o3clim
636
637 end module module_ozphys
The operational GFS currently parameterizes ozone production and destruction based on monthly mean co...