Remove older code preventing unit conversion of grid coordinates.

examples/aqm.rc

@@ -128,6 +128,11 @@ myemis_frequency: hourly
 # supported schemes are: sofiev, none
 #
 #myemis_plume_rise: sofiev
+#
+# Optionally, provide an empirical weight for the top plume layer.
+# If this option is absent, or a value < 0.0 is provided, fire
+# emissions will be distributed linearly among the plume layers
+#myemis_plume_top_fraction: 0.8
 
 # list of emission species
 # input emissions can be manipulated by adding multiple

src/aqm_cap.F90

@@ -454,20 +454,10 @@ subroutine DataInitialize(model, rc)
               rcToReturn=rc)
             return  ! bail out
           end if
-
-          call aqm_model_domain_coord_set(item, coord, de=localDe, rc=rc)
-
-          if (aqm_rc_check(rc)) then
-            call ESMF_LogSetError(ESMF_RC_INTNRL_BAD, &
-              msg="Failed to set coordinates for air quality model", &
-              line=__LINE__, &
-              file=__FILE__, &
-              rcToReturn=rc)
-            return  ! bail out
-          end if
         end do
 
       end do
+
       deallocate(minIndexPDe, maxIndexPDe, minIndexPTile, maxIndexPTile, &
         computationalLBound, computationalUBound, &
         deToTileMap, localDeToDeMap, stat=localrc)

src/model/src/PT3D_DEFN.F

@@ -55,6 +55,14 @@ FUNCTION PT3D_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME, TSTEP )
             RETURN
          END IF
 
+         ! -- merge PM maps from fire and point-source emissions
+
+         IF ( .NOT. PTMAP_INIT( ) ) THEN
+            XMSG = 'Could not merge point source mappings'
+            CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
+            SUCCESS = .FALSE.; RETURN
+         END IF
+
          ! -- initialize emission types
 
          IF ( .NOT. PT3D_FIRE_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME ) ) THEN
@@ -69,14 +77,6 @@ FUNCTION PT3D_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME, TSTEP )
             SUCCESS = .FALSE.; RETURN
          END IF
 
-         ! -- merge PM maps from fire and point-source emissions
-
-         IF ( .NOT. PTMAP_INIT( ) ) THEN
-            XMSG = 'Could not merge point source mappings'
-            CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
-            SUCCESS = .FALSE.; RETURN
-         END IF
-
          ! -- allocate emission arrays
 
          ALLOCATE ( VDEMIS_PT( NCOLS, NROWS, EMLAYS, N_SPC_PTEM ), STAT = IOS )

src/model/src/PTMAP.F

@@ -16,9 +16,6 @@ MODULE PTMAP
       INTEGER, SAVE :: N_SPC_PTEM = 0   ! merged no. of unique species
       INTEGER, SAVE :: N_SPC_PTPM = 0   ! merged no. of unique species for PM
 
-C Emission type
-c     CHARACTER( * ),  PARAMETER :: ETYPE = 'gbbepx'
-
       PRIVATE
 
       PUBLIC :: N_GSPC_EMIS, N_SPC_PTEM, N_SPC_PTPM

src/shr/aqm_emis_mod.F90

@@ -205,6 +205,7 @@ subroutine aqm_emis_src_create(config, em, rc)
       em(item) % count       = 0
       em(item) % layers      = 1
       em(item) % scalefactor = 1.0
+      em(item) % topfraction = -1.0
       em(item) % gridded     = .true.
       em(item) % sync        = .false.
       em(item) % verbose     = .false.
@@ -606,6 +607,33 @@ subroutine aqm_emis_src_init(model, em, rc)
             rcToReturn=rc)) &
             return  ! bail out
         end if
+        if (trim(em % plumerise) /= "none") then
+          call ESMF_ConfigGetAttribute(config, em % topfraction, &
+            label=trim(em % name)//"_plume_top_fraction:", default=-1.0, rc=localrc)
+          if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
+            line=__LINE__,  &
+            file=__FILE__,  &
+            rcToReturn=rc)) &
+            return  ! bail out
+          if (em % verbose) then
+            write(msgString,'(g20.8)') em % topfraction
+            call ESMF_LogWrite(trim(em % logprefix)//": "//pName &
+              //": plume_top_fraction: "//adjustl(msgString), ESMF_LOGMSG_INFO, rc=localrc)
+            if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
+              line=__LINE__,  &
+              file=__FILE__,  &
+              rcToReturn=rc)) &
+              return  ! bail out
+          end if
+          if (em % topfraction > 1.0) then
+            call ESMF_LogSetError(ESMF_RC_NOT_VALID, &
+              msg="plume top fraction must not exceed 1.0", &
+              line=__LINE__,  &
+              file=__FILE__,  &
+              rcToReturn=rc)
+            return  ! bail out
+          end if
+        end if
       case ("point-source")
         em % gridded = .false.
         ! -- get plumerise type
@@ -1621,8 +1649,8 @@ subroutine aqm_emis_grd_read(em, spcname, buffer, localDe, rc)
             return
         end select
         if (em % verbose) then
-          write(msgString, '(a12,": read",12x,": ",a16,": min/max = ",2g20.8)') &
-            em % logprefix, spcname, minval(fptr), maxval(fptr)
+          write(msgString, '(a,": read: ",a16,": min/max = ",2g20.8)') &
+            trim(em % logprefix), spcname, minval(fptr), maxval(fptr)
           call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO, rc=localrc)
           if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
             line=__LINE__,  &
@@ -1721,8 +1749,8 @@ subroutine aqm_emis_pts_read(em, spcname, buffer, ip, jp, ijmap, localDe, rc)
             em_min = min( em_min, em % rates(item) % values(n) )
             em_max = max( em_max, em % rates(item) % values(n) )
           end do
-          write(msgString, '(a12,": read",12x,": ",a16,": min/max = ",2g20.8)') &
-            em % logprefix, spcname, em_min, em_max
+          write(msgString, '(a,": read: ",a16,": min/max = ",2g20.8)') &
+            trim(em % logprefix), spcname, em_min, em_max
           call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO, rc=localrc)
           if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
             line=__LINE__,  &

src/shr/aqm_fires_mod.F90

@@ -25,6 +25,7 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc)
     integer :: c, r, l
     integer :: lev0, lev1
     real    :: Hp, pblh, th0, th1, dz
+    real    :: tfrac, w
     real(AQM_KIND_R8) :: hbl
     real(AQM_KIND_R8),    pointer :: phi(:)
     type(aqm_state_type), pointer :: state
@@ -44,16 +45,23 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc)
     ! -- get model info
     call aqm_model_get(stateIn=state, rc=localrc)
     if (aqm_rc_check(localrc, msg="Failed to retrieve model state", &
-      file=__FILE__, line=__LINE__)) return
+      file=__FILE__, line=__LINE__, rc=rc)) return
 
     ! -- get domain info
     call aqm_model_domain_get(ids=is, ide=ie, jds=js, jde=je, nl=nl, rc=localrc)
     if (aqm_rc_check(localrc, msg="Failed to retrieve grid coordinates", &
-      file=__FILE__, line=__LINE__)) return
+      file=__FILE__, line=__LINE__, rc=rc)) return
 
     nx = ie - is + 1
     ny = je - js + 1
 
+    ! -- compute layer empirical weights
+    if (aqm_rc_test((em % topfraction > 1.0), &
+      msg="Plume top fraction must be between 0.0 and 1.0", &
+      file=__FILE__, line=__LINE__, rc=rc)) return
+
+    w = min( 1.0 - em % topfraction, 1.0 )
+
     ! -- select free-troposphere vertical level
     k = 0
     do r = 1, ny
@@ -64,11 +72,11 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc)
         lev0 = minloc(phi, 1, mask=phi >= hbl)
         if (aqm_rc_test((phi(lev0) < hbl), &
           msg="Could not find first free-troposphere layer", &
-          file=__FILE__, line=__LINE__)) return
+          file=__FILE__, line=__LINE__,rc=rc)) return
         lev1 = lev0 + 1
         if (aqm_rc_test((lev1 > nl), &
           msg="Not enough vertical levels", &
-          file=__FILE__, line=__LINE__)) return
+          file=__FILE__, line=__LINE__,rc=rc)) return
 
         dz   = onebg * ( phi(lev1) - phi(lev0) )
         th0  = state % temp(c,r,lev0) * (p_ref / state % prl(c,r,lev0)) ** rcp
@@ -84,10 +92,13 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc)
 
         if (lev1 > lev0) then
           dz = 0.0
-          do l = lev0, lev1
-            profile(c,r,l) = ( phi(l) - dz ) / phi(lev1)
+          tfrac = 0.0
+          do l = lev0, lev1-1
+            profile(c,r,l) = w * ( phi(l) - dz ) / phi(lev1)
             dz = phi(l)
+            tfrac = tfrac + profile(c,r,l)
           end do
+          profile(c,r,lev1) = 1.0 - tfrac
         else
           profile(c,r,lev0) = 1.0
         end if

src/shr/aqm_internal_mod.F90

@@ -30,6 +30,7 @@ module aqm_internal_mod
     logical                        :: sync
     logical                        :: verbose
     real                           :: scalefactor
+    real                           :: topfraction
     integer(ESMF_KIND_I4)          :: layers
     integer                        :: count
     integer                        :: irec