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
51 ! Prognostic ozone.
52 integer :: nlat
53 integer :: nlev
54 integer :: ntime
55 integer :: ncf
56 real(kind_phys), allocatable :: lat(:)
57 real(kind_phys), allocatable :: pres(:)
58 real(kind_phys), allocatable :: po3(:)
59 real(kind_phys), allocatable :: time(:)
60 real(kind_phys), allocatable :: data(:,:,:,:)
61 ! Climotological ozone.
62 integer :: nlatc
63 integer :: nlevc
64 integer :: ntimec
65 real(kind_phys) :: blatc
66 real(kind_phys) :: dphiozc
67 real(kind_phys), allocatable :: pkstr(:)
68 real(kind_phys), allocatable :: pstr(:)
69 real(kind_phys), allocatable :: datac(:,:,:)
70 integer :: k1oz
71 integer :: k2oz
72 real(kind_phys) :: facoz
73 contains
74 procedure, public :: load_o3prog
75 procedure, public :: setup_o3prog
76 procedure, public :: update_o3prog
77 procedure, public :: run_o3prog_2015
78 procedure, public :: run_o3prog_2006
79 !
80 procedure, public :: load_o3clim
81 procedure, public :: update_o3clim
82 procedure, public :: run_o3clim
83 end type ty_ozphys
84
85contains
86
88 function load_o3prog(this, file, fileID) result (err_message)
89 class(ty_ozphys), intent(inout) :: this
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
96
97 ! initialize error message
98 err_message = ""
99
100 ! Get dimensions from data file
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
105 rewind(fileid)
106
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))
112
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
116
117 ! Store
118 this%pres(:) = pres4(:)
119 this%po3(:) = log(100.0*this%pres(:)) ! from mb to ln(Pa)
120 this%lat(:) = lat4(:)
121 this%time(:) = time4(:)
122 deallocate(lat4, pres4, time4)
123
124 allocate(tempin(this%nlat))
125 do i1=1,this%ntime
126 do i2=1,this%ncf
127 do i3=1,this%nlev
128 read(fileid, iostat=ierr, iomsg=err_message) tempin
129 if (ierr /= 0 ) return
130 this%data(:,i3,i2,i1) = tempin(:)
131 enddo
132 enddo
133 enddo
134 deallocate(tempin)
135 close(fileid)
136
137 end function load_o3prog
138
141 subroutine setup_o3prog(this, lat, idx1, idx2, idxh)
142 class(ty_ozphys), intent(in) :: this
143 real(kind_phys), intent(in) :: lat(:)
144 integer, intent(out) :: idx1(:), idx2(:)
145 real(kind_phys), intent(out) :: idxh(:)
146 integer :: i,j
147
148 do j=1,size(lat)
149 idx2(j) = this%nlat + 1
150 do i=1,this%nlat
151 if (lat(j) < this%lat(i)) then
152 idx2(j) = i
153 exit
154 endif
155 enddo
156 idx1(j) = max(idx2(j)-1,1)
157 idx2(j) = min(idx2(j),this%nlat)
158 if (idx2(j) .ne. idx1(j)) then
159 idxh(j) = (lat(j) - this%lat(idx1(j))) / (this%lat(idx2(j)) - this%lat(idx1(j)))
160 else
161 idxh(j) = 1.0
162 endif
163 enddo
164
165 end subroutine setup_o3prog
166
168 subroutine update_o3prog(this, idx1, idx2, idxh, rjday, idxt1, idxt2, ozpl)
169 class(ty_ozphys), intent(in) :: this
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
177
178 tx1 = (this%time(idxt2) - rjday) / (this%time(idxt2) - this%time(idxt1))
179 tx2 = 1.0 - tx1
180
181 do nc=1,this%ncf
182 do l=1,this%nlev
183 do j=1,size(ozpl(:,1,1))
184 j1 = idx1(j)
185 j2 = idx2(j)
186 tem = 1.0 - idxh(j)
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))
189 enddo
190 enddo
191 enddo
192
193 end subroutine update_o3prog
194
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)
198 class(ty_ozphys), intent(in) :: this
199 real(kind_phys), intent(in) :: &
200 con_1ovg ! Physical constant: One divided by gravitational acceleration (m-1 s2)
201 real(kind_phys), intent(in) :: &
202 dt ! Model timestep (sec)
203 real(kind_phys), intent(in), dimension(:,:) :: &
204 p, & ! Model Pressure (Pa)
205 t, & ! Model temperature (K)
206 dp ! Model layer thickness (Pa)
207 real(kind_phys), intent(in), dimension(:,:,:) :: &
208 ozpl ! Ozone forcing data
209 real(kind_phys), intent(inout), dimension(:,:) :: &
210 oz ! Ozone concentration updated by physics
211 logical, intent(in) :: do_diag
212 real(kind_phys), intent(inout), dimension(:,:), optional :: &
213 do3_dt_prd, & ! Physics tendency: production and loss effect
214 do3_dt_ozmx, & ! Physics tendency: ozone mixing ratio effect
215 do3_dt_temp, & ! Physics tendency: temperature effect
216 do3_dt_ohoz ! Physics tendency: overhead ozone effect
217
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
225
226 ! Dimensions
227 ncol = size(p,1)
228 nlev = size(p,2)
229
230 ! Temporaries
231 ozi = oz
232
233 colo3(:,nlev+1) = 0.0
234 coloz(:,nlev+1) = 0.0
235
236 do ilev=nlev,1,-1
237 pmin = 1.0e10
238 pmax = -1.0e10
239
240 do icol=1,ncol
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
245 enddo
246 kmax = 1
247 kmin = 1
248 do k=1,this%nlev-1
249 if (pmin < this%po3(k)) kmax = k
250 if (pmax < this%po3(k)) kmin = k
251 enddo
252 !
253 do k=kmin,kmax
254 temp = 1.0 / (this%po3(k) - this%po3(k+1))
255 do icol=1,ncol
256 flg(icol) = .false.
257 if (wk1(icol) < this%po3(k) .and. wk1(icol) >= this%po3(k+1)) then
258 flg(icol) = .true.
259 wk2(icol) = (wk1(icol) - this%po3(k+1)) * temp
260 wk3(icol) = 1.0 - wk2(icol)
261 endif
262 enddo
263 do icf=1,this%ncf
264 do icol=1,ncol
265 if (flg(icol)) then
266 prod(icol,icf) = wk2(icol) * ozpl(icol,k,icf) + wk3(icol) * ozpl(icol,k+1,icf)
267 endif
268 enddo
269 enddo
270 enddo
271
272 do icf=1,this%ncf
273 do icol=1,ncol
274 if (wk1(icol) < this%po3(this%nlev)) then
275 prod(icol,icf) = ozpl(icol,this%nlev,icf)
276 endif
277 if (wk1(icol) >= this%po3(1)) then
278 prod(icol,icf) = ozpl(icol,1,icf)
279 endif
280 enddo
281 enddo
282 do icol=1,ncol
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)
286 enddo
287 do icol=1,ncol
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)
293 enddo
294
295 ! Diagnostics (optional)
296 if (do_diag) then
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
301 endif
302 enddo
303
304 return
305 end subroutine run_o3prog_2015
306
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)
310 class(ty_ozphys), intent(in) :: this
311 real(kind_phys), intent(in) :: &
312 con_1ovg ! Physical constant: One divided by gravitational acceleration (m-1 s2)
313 real(kind_phys), intent(in) :: &
314 dt ! Model timestep (sec)
315 real(kind_phys), intent(in), dimension(:,:) :: &
316 p, & ! Model Pressure (Pa)
317 t, & ! Model temperature (K)
318 dp ! Model layer thickness (Pa)
319 real(kind_phys), intent(in), dimension(:,:,:) :: &
320 ozpl ! Ozone forcing data
321 real(kind_phys), intent(inout), dimension(:,:) :: &
322 oz ! Ozone concentration updated by physics
323 logical, intent(in) :: do_diag
324 real(kind_phys), intent(inout), dimension(:,:) :: &
325 do3_dt_prd, & ! Physics tendency: production and loss effect
326 do3_dt_ozmx, & ! Physics tendency: ozone mixing ratio effect
327 do3_dt_temp, & ! Physics tendency: temperature effect
328 do3_dt_ohoz ! Physics tendency: overhead ozone effect
329
330 ! Locals
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
338
339 ! Dimensions
340 ncol = size(p,1)
341 nlev = size(p,2)
342
343 ! Temporaries
344 ozi = oz
345
347 if (this%ncf > 2) then
348 colo3(:,nlev+1) = 0.0
349 do ilev=nlev,1,-1
350 do icol=1,ncol
351 colo3(icol,ilev) = colo3(icol,ilev+1) + ozi(icol,ilev) * dp(icol,ilev) * con_1ovg
352 enddo
353 enddo
354 endif
355
357 do ilev=1,nlev
358 pmin = 1.0e10
359 pmax = -1.0e10
360
361 do icol=1,ncol
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
366 enddo
367 kmax = 1
368 kmin = 1
369 do k=1,this%nlev-1
370 if (pmin < this%po3(k)) kmax = k
371 if (pmax < this%po3(k)) kmin = k
372 enddo
373
374 do k=kmin,kmax
375 temp = 1.0 / (this%po3(k) - this%po3(k+1))
376 do icol=1,ncol
377 flg(icol) = .false.
378 if (wk1(icol) < this%po3(k) .and. wk1(icol) >= this%po3(k+1)) then
379 flg(icol) = .true.
380 wk2(icol) = (wk1(icol) - this%po3(k+1)) * temp
381 wk3(icol) = 1.0 - wk2(icol)
382 endif
383 enddo
384 do icf=1,this%ncf
385 do icol=1,ncol
386 if (flg(icol)) then
387 prod(icol,icf) = wk2(icol) * ozpl(icol,k,icf) + wk3(icol) * ozpl(icol,k+1,icf)
388 endif
389 enddo
390 enddo
391 enddo
392
393 do icf=1,this%ncf
394 do icol=1,ncol
395 if (wk1(icol) < this%po3(this%nlev)) then
396 prod(icol,icf) = ozpl(icol,this%nlev,icf)
397 endif
398 if (wk1(icol) >= this%po3(1)) then
399 prod(icol,icf) = ozpl(icol,1,icf)
400 endif
401 enddo
402 enddo
403
404 if (this%ncf == 2) then
405 do icol=1,ncol
406 ozib(icol) = ozi(icol,ilev)
407 oz(icol,ilev) = (ozib(icol) + prod(icol,1)*dt) / (1.0 + prod(icol,2)*dt)
408 enddo
409 endif
410
411 if (this%ncf == 4) then
412 do icol=1,ncol
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)
416 enddo
417 endif
418
419 ! Diagnostics (optional)
420 if (do_diag) then
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
425 endif
426 enddo
427
428 return
429 end subroutine run_o3prog_2006
430
432 subroutine run_o3clim(this, lat, prslk, con_pi, oz)
433 class(ty_ozphys), intent(in) :: this
434 real(kind_phys), intent(in) :: &
435 con_pi ! Physics constant: Pi
436 real(kind_phys), intent(in), dimension(:) :: &
437 lat ! Grid latitude
438 real(kind_phys), intent(in), dimension(:,:) :: &
439 prslk ! Exner function
440 real(kind_phys), intent(out), dimension(:,:) :: &
441 oz ! Ozone concentration updated by climotology
442
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(:)
446 logical :: top_at_1
447
448 ncol = size(prslk(:,1))
449 nlev = size(prslk(1,:))
450 allocate(o3i(ncol, this%nlevc),wk1(ncol))
451
452 ! What is vertical ordering?
453 top_at_1 = (prslk(1,1) .lt. prslk(1, nlev))
454
455 elte = this%blatc + (this%nlatc-1)*this%dphiozc
456
457 do icol = 1, ncol
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
461 j1 = tem1
462 j2 = j1 + 1
463 tem1 = tem1 - j1
464 elseif (deglat <= this%blatc) then
465 j1 = 1
466 j2 = 1
467 tem1 = 1.0
468 elseif (deglat >= elte) then
469 j1 = this%nlatc
470 j2 = this%nlatc
471 tem1 = 1.0
472 endif
473
474 tem2 = 1.0 - tem1
475 do j = 1, this%nlevc
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)
479 enddo
480 enddo
481
482 do ilev = 1, nlev
483 ll = ilev
484 if (.not. top_at_1) ll = nlev - ilev + 1
485
486 do icol = 1, ncol
487 wk1(icol) = prslk(icol,ll)
488 enddo
489
490 do j = 1, this%nlevc-1
491 temp = 1.0 / (this%pkstr(j+1) - this%pkstr(j))
492 do icol = 1, ncol
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)
496 endif
497 enddo
498 enddo
499
500 do icol = 1, ncol
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)
503 enddo
504 enddo
505
506 return
507 end subroutine run_o3clim
508
510 function load_o3clim(this, file, fileID) result (err_message)
511 class(ty_ozphys), intent(inout) :: this
512 integer, intent(in) :: fileid
513 character(len=*), intent(in) :: file
514 character(len=128) :: err_message
515
516 ! Locals
517 real(kind=4) :: blatc4
518 integer :: ilev, ilat, imo, ierr
519 real(kind=4), allocatable :: o3clim4(:,:,:), pstr4(:)
520 integer, allocatable :: imond(:), ilatt(:,:)
521
522 ! initialize error message
523 err_message = ""
524
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
529
530101 if (this%nlevc < 10 .or. this%nlevc > 100) then
531 rewind(fileid)
532 this%nlevc = 17
533 this%nlatc = 18
534 this%blatc = -85.0
535 else
536 this%blatc = blatc4
537 endif
538 this%nlat = 2
539 this%nlev = 1
540 this%ntimec = 1
541 this%ncf = 0
542 this%dphiozc = -(this%blatc+this%blatc)/(this%nlatc-1)
543
544 allocate (o3clim4(this%nlatc,this%nlevc,12), pstr4(this%nlevc), imond(12), ilatt(this%nlatc,12) )
545
546 allocate (this%pkstr(this%nlevc), this%pstr(this%nlevc), this%datac(this%nlatc,this%nlevc,12))
547 if ( this%nlevc == 17 ) then ! For the operational ozone climatology
548 do ilev = 1, this%nlevc
549 read (fileid,15,iostat=ierr,iomsg=err_message) pstr4(ilev)
550 if (ierr /= 0 ) return
55115 format(f10.3)
552 enddo
553
554 do imo = 1, 12
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
56120 format(6x,10f6.2)
562 enddo
563 enddo
564 else ! For newer ozone climatology
565 read (fileid)
566 do ilev = 1, this%nlevc
567 read (fileid) pstr4(ilev)
568 enddo
569
570 do imo = 1, 12
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
574 enddo
575 enddo
576 endif ! end if_this%nlevc_block
577
578 do imo = 1, 12
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
582 enddo
583 enddo
584 enddo
585
586 do ilev = 1, this%nlevc
587 this%pstr(ilev) = pstr4(ilev)
588 this%pkstr(ilev) = fpkapx(this%pstr(ilev)*100.0)
589 enddo
590
591 end function load_o3clim
592
595 subroutine update_o3clim(this, imon, iday, ihour, loz1st)
596 class(ty_ozphys), intent(inout) :: this
597 integer, intent(in) :: imon, iday, ihour
598 logical, intent(in) :: loz1st
599
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/)
602 logical :: change
603
604 midmon = mdays(imon)/2 + 1
605 change = loz1st .or. ( (iday==midmon) .and. (ihour==0) )
606
607 if ( change ) then
608 if ( iday < midmon ) then
609 this%k1oz = mod(imon+10, 12) + 1
610 midm = mdays(this%k1oz)/2 + 1
611 this%k2oz = imon
612 midp = mdays(this%k1oz) + midmon
613 else
614 this%k1oz = imon
615 midm = midmon
616 this%k2oz = mod(imon, 12) + 1
617 midp = mdays(this%k2oz)/2 + 1 + mdays(this%k1oz)
618 endif
619 endif
620
621 if (iday < midmon) then
622 id = iday + mdays(this%k1oz)
623 else
624 id = iday
625 endif
626
627 this%facoz = float(id - midm) / float(midp - midm)
628
629 end subroutine update_o3clim
630
631 end module module_ozphys
The operational GFS currently parameterizes ozone production and destruction based on monthly mean co...
Derived type containing data and procedures needed by ozone photochemistry parameterization Note All ...