Reorganize ozone physics into CCPP phases

physics/GFS_phys_time_vary.fv3.F90

@@ -2,20 +2,19 @@
 !!  Contains code related to GFS physics suite setup (physics part of time_vary_step)
 
 !>\defgroup mod_GFS_phys_time_vary GFS Physics Time Update
-!! This module contains GFS physics time vary subroutines including ozone, stratospheric water vapor,
+!! This module contains GFS physics time vary subroutines including stratospheric water vapor,
 !! aerosol, IN&CCN and surface properties updates.
    module GFS_phys_time_vary
 
 #ifdef _OPENMP
       use omp_lib
 #endif
 
-      use machine, only : kind_phys
+      use machine, only : kind_phys, kind_dbl_prec, kind_sngl_prec
 
       use mersenne_twister, only: random_setseed, random_number
 
-      use ozne_def, only : levozp, oz_coeff, oz_lat, oz_pres, oz_time, ozplin
-      use ozinterp, only : read_o3data, setindxoz, ozinterpol
+      use module_ozphys, only: ty_ozphys
 
       use h2o_def,   only : levh2o, h2o_coeff, h2o_lat, h2o_pres, h2o_time, h2oplin
       use h2ointerp, only : read_h2odata, setindxh2o, h2ointerpol
@@ -85,7 +84,7 @@ end subroutine copy_error
       subroutine GFS_phys_time_vary_init (                                                         &
               me, master, ntoz, h2o_phys, iaerclm, iccn, iaermdl, iflip, im, levs,                 &
               nx, ny, idate, xlat_d, xlon_d,                                                       &
-              jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl,fhour,          &
+              jindx1_o3, jindx2_o3, ddy_o3, jindx1_h, jindx2_h, ddy_h, h2opl,fhour,                &
               jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm,            &
               jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, imap, jmap,              &
               do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau,                            &
@@ -98,7 +97,7 @@ subroutine GFS_phys_time_vary_init (
               smcwtdxy, deeprechxy, rechxy, snowxy, snicexy, snliqxy, tsnoxy , smoiseq, zsnsoxy,   &
               slc, smc, stc, tsfcl, snowd, canopy, tg3, stype, con_t0c, lsm_cold_start, nthrds,    &
               lkm, use_lake_model, lakefrac, lakedepth, iopt_lake, iopt_lake_clm, iopt_lake_flake, &
-              lakefrac_threshold, lakedepth_threshold, errmsg, errflg)
+              lakefrac_threshold, lakedepth_threshold, ozphys, errmsg, errflg)
 
          implicit none
 
@@ -115,7 +114,8 @@ subroutine GFS_phys_time_vary_init (
 
          integer,              intent(inout) :: jindx1_o3(:), jindx2_o3(:), jindx1_h(:), jindx2_h(:)
          real(kind_phys),      intent(inout) :: ddy_o3(:),  ddy_h(:)
-         real(kind_phys),      intent(in)    :: ozpl(:,:,:), h2opl(:,:,:)
+         real(kind_phys),      intent(in)    :: h2opl(:,:,:)
+
          integer,              intent(inout) :: jindx1_aer(:), jindx2_aer(:), iindx1_aer(:), iindx2_aer(:)
          real(kind_phys),      intent(inout) :: ddy_aer(:), ddx_aer(:)
          real(kind_phys),      intent(out)   :: aer_nm(:,:,:)
@@ -132,6 +132,7 @@ subroutine GFS_phys_time_vary_init (
          real(kind_phys),      intent(in)    :: min_seaice, fice(:)
          real(kind_phys),      intent(in)    :: landfrac(:)
          real(kind_phys),      intent(inout) :: weasd(:)
+         type(ty_ozphys),      intent(in)    :: ozphys
 
          ! NoahMP - only allocated when NoahMP is used
          integer, intent(in) :: lsoil, lsnow_lsm_lbound, lsnow_lsm_ubound
@@ -224,7 +225,6 @@ subroutine GFS_phys_time_vary_init (
 !$OMP parallel num_threads(nthrds) default(none)                                    &
 !$OMP          shared (me,master,ntoz,h2o_phys,im,nx,ny,levs,idate)                 &
 !$OMP          shared (xlat_d,xlon_d,imap,jmap,errmsg,errflg)                       &
-!$OMP          shared (levozp,oz_coeff,oz_pres,ozpl)                                &
 !$OMP          shared (levh2o,h2o_coeff,h2o_pres,h2opl)                             &
 !$OMP          shared (iamin, iamax, jamin, jamax, lsm_noahmp)                      &
 !$OMP          shared (iaerclm,iaermdl,ntrcaer,aer_nm,iflip,iccn)                   &
@@ -234,41 +234,18 @@ subroutine GFS_phys_time_vary_init (
 !$OMP          shared (do_ugwp_v1,jindx1_tau,jindx2_tau,ddy_j1tau,ddy_j2tau)        &
 !$OMP          shared (isot,ivegsrc,nlunit,sncovr,sncovr_ice,lsm,lsm_ruc)           &
 !$OMP          shared (min_seaice,fice,landfrac,vtype,weasd,snupx,salp_data)        &
+!$OMP          shared (ozphys)                                                      &
 !$OMP          private (ix,i,j,rsnow,vegtyp,myerrmsg,myerrflg)
 
 !$OMP sections
 
-!$OMP section
-!> - Call read_o3data() to read ozone data
-       need_o3data: if(ntoz > 0) then
-         call read_o3data (ntoz, me, master)
-
-         ! Consistency check that the hardcoded values for levozp and
-         ! oz_coeff in GFS_typedefs.F90 match what is set by read_o3data
-         ! in GFS_typedefs.F90: allocate (Tbd%ozpl (IM,levozp,oz_coeff))
-         if (size(ozpl, dim=2).ne.levozp) then
-            myerrflg = 1
-            write(myerrmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ",    &
-                  "levozp from read_o3data does not match value in GFS_typedefs.F90: ", &
-                  levozp, " /= ", size(ozpl, dim=2)
-            call copy_error(myerrmsg, myerrflg, errmsg, errflg)
-         end if
-         if (size(ozpl, dim=3).ne.oz_coeff) then
-            myerrflg = 1
-            write(myerrmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ",      &
-                  "oz_coeff from read_o3data does not match value in GFS_typedefs.F90: ", &
-                  oz_coeff, " /= ", size(ozpl, dim=3)
-            call copy_error(myerrmsg, myerrflg, errmsg, errflg)
-         end if
-       endif need_o3data
-
 !$OMP section
 !> - Call read_h2odata() to read stratospheric water vapor data
        need_h2odata: if(h2o_phys) then
          call read_h2odata (h2o_phys, me, master)
 
          ! Consistency check that the hardcoded values for levh2o and
-         ! h2o_coeff in GFS_typedefs.F90 match what is set by read_o3data
+         ! h2o_coeff in GFS_typedefs.F90 match what is set by read_h2odata
          ! in GFS_typedefs.F90: allocate (Tbd%h2opl (IM,levh2o,h2o_coeff))
          if (size(h2opl, dim=2).ne.levh2o) then
             write(myerrmsg,'(2a,i0,a,i0)') "Value error in GFS_phys_time_vary_init: ",     &
@@ -348,9 +325,9 @@ subroutine GFS_phys_time_vary_init (
 !$OMP sections
 
 !$OMP section
-!> - Call setindxoz() to initialize ozone data
+!> - Setup spatial interpolation indices for ozone physics.
          if (ntoz > 0) then
-           call setindxoz (im, xlat_d, jindx1_o3, jindx2_o3, ddy_o3)
+           call ozphys%setup_o3prog(xlat_d, jindx1_o3, jindx2_o3, ddy_o3)
          endif
 
 !$OMP section
@@ -794,7 +771,7 @@ subroutine GFS_phys_time_vary_timestep_init (
             lakefrac, min_seaice, min_lakeice, smc, slc, stc, smois, sh2o, tslb, tiice, tg3, tref,  &
             tsfc, tsfco, tisfc, hice, fice, facsf, facwf, alvsf, alvwf, alnsf, alnwf, zorli, zorll, &
             zorlo, weasd, slope, snoalb, canopy, vfrac, vtype, stype,scolor, shdmin, shdmax, snowd, &
-            cv, cvb, cvt, oro, oro_uf, xlat_d, xlon_d, slmsk, landfrac,                             &
+            cv, cvb, cvt, oro, oro_uf, xlat_d, xlon_d, slmsk, landfrac, ozphys,                     &
             do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau, tau_amf, errmsg, errflg)
 
          implicit none
@@ -824,6 +801,7 @@ subroutine GFS_phys_time_vary_timestep_init (
          integer,              intent(in)    :: jindx1_tau(:), jindx2_tau(:)
          real(kind_phys),      intent(in)    :: ddy_j1tau(:), ddy_j2tau(:)
          real(kind_phys),      intent(inout) :: tau_amf(:)
+         type(ty_ozphys),      intent(in)    :: ozphys
 
          ! For gcycle only
          integer,              intent(in)    :: nthrds, nx, ny, nsst, tile_num, nlunit, lsoil
@@ -846,10 +824,13 @@ subroutine GFS_phys_time_vary_timestep_init (
          integer,              intent(out)   :: errflg
 
          ! Local variables
-         integer :: i, j, k, iseed, iskip, ix
-         real(kind=kind_phys) :: wrk(1)
-         real(kind=kind_phys) :: rannie(cny)
-         real(kind=kind_phys) :: rndval(cnx*cny*nrcm)
+         integer :: i, j, k, iseed, iskip, ix, idat(8), jdat(8), iday, j1, j2, nc, n1, n2, jdow,     &
+              jdoy, jday, w3kindreal, w3kindint
+         real(kind_phys) :: wrk(1), tem, tx1, tx2, rjday
+         real(kind_phys) :: rannie(cny)
+         real(kind_phys) :: rndval(cnx*cny*nrcm)
+         real(kind_dbl_prec)  :: rinc(5)
+         real(kind_sngl_prec) :: rinc4(5)
 
          ! Initialize CCPP error handling variables
          errmsg = ''
@@ -869,7 +850,8 @@ subroutine GFS_phys_time_vary_timestep_init (
 !$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)                                           &
+!$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
@@ -920,11 +902,41 @@ subroutine GFS_phys_time_vary_timestep_init (
          endif  ! imfdeepcnv, cal_re, random_clds
 
 !$OMP section
-!> - Call ozinterpol() to make ozone interpolation
+         !> - 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 ozinterpol (me, im, idate, fhour, &
-                            jindx1_o3, jindx2_o3, &
-                            ozpl, ddy_o3)
+            call ozphys%update_o3prog(jindx1_o3, jindx2_o3, ddy_o3, rjday, n1, n2, ozpl)
          endif
 
 !$OMP section
@@ -1024,12 +1036,6 @@ subroutine GFS_phys_time_vary_finalize(errmsg, errflg)
 
          if (.not.is_initialized) return
 
-         ! Deallocate ozone arrays
-         if (allocated(oz_lat)  ) deallocate(oz_lat)
-         if (allocated(oz_pres) ) deallocate(oz_pres)
-         if (allocated(oz_time) ) deallocate(oz_time)
-         if (allocated(ozplin)  ) deallocate(ozplin)
-
          ! Deallocate h2o arrays
          if (allocated(h2o_lat) ) deallocate(h2o_lat)
          if (allocated(h2o_pres)) deallocate(h2o_pres)

physics/GFS_phys_time_vary.fv3.meta

@@ -2,7 +2,7 @@
   name = GFS_phys_time_vary
   type = scheme
   dependencies = aerclm_def.F,aerinterp.F90,gcycle.F90,h2o_def.f,h2ointerp.f90,iccn_def.F,iccninterp.F90,machine.F,mersenne_twister.f
-  dependencies = namelist_soilveg.f,set_soilveg.f,ozinterp.f90,ozne_def.f,sfcsub.F,cires_tauamf_data.F90,noahmp_tables.f90
+  dependencies = namelist_soilveg.f,set_soilveg.f,sfcsub.F,cires_tauamf_data.F90,noahmp_tables.f90,module_ozphys.F90
 
 ########################################################################
 [ccpp-arg-table]
@@ -138,14 +138,6 @@
   type = real
   kind = kind_phys
   intent = inout
-[ozpl]
-  standard_name = ozone_forcing
-  long_name = ozone forcing data
-  units = mixed
-  dimensions = (horizontal_dimension,vertical_dimension_of_ozone_forcing_data,number_of_coefficients_in_ozone_forcing_data)
-  type = real
-  kind = kind_phys
-  intent = in
 [jindx1_h]
   standard_name = lower_latitude_index_of_stratospheric_water_vapor_forcing_for_interpolation
   long_name = interpolation low index for stratospheric water vapor
@@ -969,6 +961,13 @@
   type = real
   kind = kind_phys
   intent = in
+[ozphys]
+  standard_name = dataset_for_ozone_physics
+  long_name = dataset for NRL ozone physics
+  units = mixed
+  dimensions = ()
+  type = ty_ozphys
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP
@@ -1206,7 +1205,7 @@
   standard_name = ozone_forcing
   long_name = ozone forcing data
   units = mixed
-  dimensions = (horizontal_dimension,vertical_dimension_of_ozone_forcing_data,number_of_coefficients_in_ozone_forcing_data)
+  dimensions = (horizontal_dimension,vertical_dimension_of_ozone_forcing_data,number_of_coefficients_in_ozone_data)
   type = real
   kind = kind_phys
   intent = inout
@@ -1942,6 +1941,13 @@
   type = real
   kind = kind_phys
   intent = inout
+[ozphys]
+  standard_name = dataset_for_ozone_physics
+  long_name = dataset for NRL ozone physics
+  units = mixed
+  dimensions = ()
+  type = ty_ozphys
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP