Skip to content

Commit

Permalink
Some polishing. Merge 2006 ozone into module_ozphys
Browse files Browse the repository at this point in the history
  • Loading branch information
dustinswales committed Sep 28, 2023
1 parent d0a4bfd commit 385ef4e
Show file tree
Hide file tree
Showing 5 changed files with 158 additions and 460 deletions.
69 changes: 35 additions & 34 deletions physics/GFS_phys_time_vary.fv3.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -802,47 +802,16 @@ subroutine GFS_phys_time_vary_timestep_init (
return
end if

!> - Compute temporal interpolation indices for updating gas concentrations.
idat=0
idat(1)=idate(4)
idat(2)=idate(2)
idat(3)=idate(3)
idat(5)=idate(1)
rinc=0.
rinc(2)=fhour
call w3kind(w3kindreal,w3kindint)
if(w3kindreal==4) then
rinc4=rinc
CALL w3movdat(rinc4,idat,jdat)
else
CALL w3movdat(rinc,idat,jdat)
endif
jdow = 0
jdoy = 0
jday = 0
call w3doxdat(jdat,jdow,jdoy,jday)
rjday = jdoy + jdat(5) / 24.
if (rjday < ozphys%time(1)) rjday = rjday + 365.

n2 = ozphys%ntime + 1
do j=2,ozphys%ntime
if (rjday < ozphys%time(j)) then
n2 = j
exit
endif
enddo
n1 = n2 - 1
if (n2 > ozphys%ntime) n2 = n2 - ozphys%ntime

!$OMP parallel num_threads(nthrds) default(none) &
!$OMP shared(kdt,nsswr,lsswr,clstp,imfdeepcnv,cal_pre,random_clds) &
!$OMP shared(fhswr,fhour,seed0,cnx,cny,nrcm,wrk,rannie,rndval) &
!$OMP shared(rann,im,isc,jsc,imap,jmap,ntoz,me,idate,jindx1_o3,jindx2_o3) &
!$OMP shared(ozpl,ddy_o3,h2o_phys,jindx1_h,jindx2_h,h2opl,ddy_h,iaerclm,master) &
!$OMP shared(levs,prsl,iccn,jindx1_ci,jindx2_ci,ddy_ci,iindx1_ci,iindx2_ci) &
!$OMP shared(ddx_ci,in_nm,ccn_nm,do_ugwp_v1,jindx1_tau,jindx2_tau,ddy_j1tau) &
!$OMP shared(ddy_j2tau,tau_amf,iflip,ozphys) &
!$OMP private(iseed,iskip,i,j,k,rjday,n1,n2)
!$OMP shared(ddy_j2tau,tau_amf,iflip,ozphys,rjday,n1,n2,idat,jdat,rinc,rinc4) &
!$OMP shared(w3kindreal,w3kindint,jdow,jdoy,jday) &
!$OMP private(iseed,iskip,i,j,k)

!$OMP sections

Expand Down Expand Up @@ -892,6 +861,38 @@ subroutine GFS_phys_time_vary_timestep_init (
endif ! imfdeepcnv, cal_re, random_clds

!$OMP section
!> - Compute temporal interpolation indices for updating gas concentrations.
idat=0
idat(1)=idate(4)
idat(2)=idate(2)
idat(3)=idate(3)
idat(5)=idate(1)
rinc=0.
rinc(2)=fhour
call w3kind(w3kindreal,w3kindint)
if(w3kindreal==4) then
rinc4=rinc
CALL w3movdat(rinc4,idat,jdat)
else
CALL w3movdat(rinc,idat,jdat)
endif
jdow = 0
jdoy = 0
jday = 0
call w3doxdat(jdat,jdow,jdoy,jday)
rjday = jdoy + jdat(5) / 24.
if (rjday < ozphys%time(1)) rjday = rjday + 365.

n2 = ozphys%ntime + 1
do j=2,ozphys%ntime
if (rjday < ozphys%time(j)) then
n2 = j
exit
endif
enddo
n1 = n2 - 1
if (n2 > ozphys%ntime) n2 = n2 - ozphys%ntime

!> - Update ozone concentration.
if (ntoz > 0) then
call ozphys%update_o3prog(jindx1_o3, jindx2_o3, ddy_o3, rjday, n1, n2, ozpl)
Expand Down
3 changes: 2 additions & 1 deletion physics/GFS_suite_stateout_update.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,8 @@ subroutine GFS_suite_stateout_update_run (im, levs, ntrac, dtp, tgrs, ugrs, vgrs
do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz)
endif
if (oz_phys_2006) then
call ozphys%run_o3prog_2006()
call ozphys%run_o3prog_2006(con_1ovg, dtp, prsl, gt0, dp, ozpl, oz0, do3_dt_prd, &
do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz)
endif

! If using Ferrier-Aligo microphysics, set bounds on the mass-weighted rime factor.
Expand Down
127 changes: 121 additions & 6 deletions physics/module_ozphys.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,12 +162,12 @@ end subroutine update_o3prog
! #########################################################################################
! Procedure (type-bound) for NRL prognostic ozone (2015).
! #########################################################################################
subroutine run_o3prog_2015(this, con_1ovg, dt, p, t, dp, ozpl, oz, do3_dt_prd, do3_dt_ozmx, &
do3_dt_temp, do3_dt_ohoz)
subroutine run_o3prog_2015(this, con_1ovg, dt, p, t, dp, ozpl, oz, do3_dt_prd, &
do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz)
class(ty_ozphys), intent(in) :: this
real(kind_phys),intent(in) :: &
real(kind_phys), intent(in) :: &
con_1ovg ! Physical constant: One divided by gravitational acceleration (m-1 s2)
real(kind_phys), intent(in) :: &
real(kind_phys), intent(in) :: &
dt ! Model timestep (sec)
real(kind_phys), intent(in), dimension(:,:) :: &
p, & ! Model Pressure (Pa)
Expand Down Expand Up @@ -253,7 +253,7 @@ subroutine run_o3prog_2015(this, con_1ovg, dt, p, t, dp, ozpl, oz, do3_dt_prd, d
prod(iCol,2) = min(prod(iCol,2), 0.0)
enddo
do iCol=1,nCol
ozib(iCol) = ozi(iCol,iLev) ! no filling
ozib(iCol) = ozi(iCol,iLev)
tem = prod(iCol,1) - prod(iCol,2) * prod(iCol,6) &
+ prod(iCol,3) * (t(iCol,iLev) - prod(iCol,5)) &
+ prod(iCol,4) * (colo3(iCol,iLev)-coloz(iCol,iLev))
Expand All @@ -273,8 +273,123 @@ end subroutine run_o3prog_2015
! #########################################################################################
! Procedure (type-bound) for NRL prognostic ozone (2006).
! #########################################################################################
subroutine run_o3prog_2006(this)
subroutine run_o3prog_2006(this, con_1ovg, dt, p, t, dp, ozpl, oz, do3_dt_prd, &
do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz)
class(ty_ozphys), intent(in) :: this
real(kind_phys), intent(in) :: &
con_1ovg ! Physical constant: One divided by gravitational acceleration (m-1 s2)
real(kind_phys), intent(in) :: &
dt ! Model timestep (sec)
real(kind_phys), intent(in), dimension(:,:) :: &
p, & ! Model Pressure (Pa)
t, & ! Model temperature (K)
dp ! Model layer thickness (Pa)
real(kind_phys), intent(in), dimension(:,:,:) :: &
ozpl ! Ozone forcing data
real(kind_phys), intent(inout), dimension(:,:) :: &
oz ! Ozone concentration updated by physics
real(kind_phys), intent(inout), dimension(:,:), pointer, optional :: &
do3_dt_prd, & ! Physics tendency: production and loss effect
do3_dt_ozmx, & ! Physics tendency: ozone mixing ratio effect
do3_dt_temp, & ! Physics tendency: temperature effect
do3_dt_ohoz ! Physics tendency: overhead ozone effect

! Locals
integer :: k, kmax, kmin, iLev, iCol, nCol, nLev, iCf
logical, dimension(size(p,1)) :: flg
real(kind_phys) :: pmax, pmin, tem, temp
real(kind_phys), dimension(size(p,1)) :: wk1, wk2, wk3, ozib
real(kind_phys), dimension(size(p,1),this%ncf) :: prod
real(kind_phys), dimension(size(p,1),size(p,2)) :: ozi
real(kind_phys), dimension(size(p,1),size(p,2)+1) :: colo3, coloz

! Dimensions
nCol = size(p,1)
nLev = size(p,2)

! Temporaries
ozi = oz

!> - Calculate vertical integrated column ozone values.
if (this%ncf > 2) then
colo3(:,nLev+1) = 0.0
do iLev=nLev,1,-1
do iCol=1,nCol
colo3(iCol,iLev) = colo3(iCol,iLev+1) + ozi(iCol,iLev) * dp(iCol,iLev) * con_1ovg
enddo
enddo
endif

!> - Apply vertically linear interpolation to the ozone coefficients.
do iLev=1,nLev
pmin = 1.0e10
pmax = -1.0e10

do iCol=1,nCol
wk1(iCol) = log(p(iCol,iLev))
pmin = min(wk1(iCol), pmin)
pmax = max(wk1(iCol), pmax)
prod(iCol,:) = 0._kind_phys
enddo
kmax = 1
kmin = 1
do k=1,this%nlev-1
if (pmin < this%po3(k)) kmax = k
if (pmax < this%po3(k)) kmin = k
enddo

do k=kmin,kmax
temp = 1.0 / (this%po3(k) - this%po3(k+1))
do iCol=1,nCol
flg(iCol) = .false.
if (wk1(iCol) < this%po3(k) .and. wk1(iCol) >= this%po3(k+1)) then
flg(iCol) = .true.
wk2(iCol) = (wk1(iCol) - this%po3(k+1)) * temp
wk3(iCol) = 1.0 - wk2(iCol)
endif
enddo
do iCf=1,this%ncf
do iCol=1,nCol
if (flg(iCol)) then
prod(iCol,iCf) = wk2(iCol) * ozpl(iCol,k,iCf) + wk3(iCol) * ozpl(iCol,k+1,iCf)
endif
enddo
enddo
enddo

do iCf=1,this%ncf
do iCol=1,nCol
if (wk1(iCol) < this%po3(this%nlev)) then
prod(iCol,iCf) = ozpl(iCol,this%nlev,iCf)
endif
if (wk1(iCol) >= this%po3(1)) then
prod(iCol,iCf) = ozpl(iCol,1,iCf)
endif
enddo
enddo

if (this%ncf == 2) then
do iCol=1,nCol
ozib(iCol) = ozi(iCol,iLev)
oz(iCol,iLev) = (ozib(iCol) + prod(iCol,1)*dt) / (1.0 + prod(iCol,2)*dt)
enddo
endif

if (this%ncf == 4) then
do iCol=1,nCol
ozib(iCol) = ozi(iCol,iLev)
tem = prod(iCol,1) + prod(iCol,3)*t(iCol,iLev) + prod(iCol,4)*colo3(iCol,iLev+1)
oz(iCol,iLev) = (ozib(iCol) + tem*dt) / (1.0 + prod(iCol,2)*dt)
enddo
endif
! Diagnostics (optional)
if (associated(do3_dt_prd)) do3_dt_prd(:,iLev) = prod(:,1)*dt
if (associated(do3_dt_ozmx)) do3_dt_ozmx(:,iLev) = (oz(:,iLev) - ozib(:))
if (associated(do3_dt_temp)) do3_dt_temp(:,iLev) = prod(:,3) * t(:,iLev) * dt
if (associated(do3_dt_ohoz)) do3_dt_ohoz(:,iLev) = prod(:,4) * colo3(:,iLev) * dt

enddo

return
end subroutine run_o3prog_2006

Expand Down
Loading

0 comments on commit 385ef4e

Please sign in to comment.