Merge pull request #131 from GEOS-ESM/hotfix/mathomp4/#130-fix-klid-i…

…ssue

Fixes #130. Initialize allocatable variables in Process Library

CHANGELOG.md

@@ -7,12 +7,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [2.0.6] - 2021-04-28
+
+### Fixed
+
+- Initialize allocatable variables in Process Library. Fixes #130
+
 ## [2.0.5] - 2021-03-14
 
 ### Added
 
 - Added AMIP.20C ExtData configs to allow AMIP GOCART runs to work before Y2000 (during the transition from HFED to QFED)
-  - Note 1: This is not a *new* scenario but rather a stopgap until Extdata is updated to allow time ranges to be specified. 
+  - Note 1: This is not a *new* scenario but rather a stopgap until Extdata is updated to allow time ranges to be specified.
   - Note 2: Temporarily, this will allow runs before Y2000 using magic-data scripting in `gcm_run.j` (a la MERRA2 in GOCART1).
 
 ### Fixed

CMakeLists.txt

@@ -4,7 +4,7 @@ cmake_policy (SET CMP0054 NEW)
 
 project (
         GOCART
-        VERSION 2.0.5
+        VERSION 2.0.6
         LANGUAGES Fortran CXX C)  # Note - CXX is required for ESMF
 
 if ("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}")

Process_Library/GOCART2G_Process.F90

@@ -661,8 +661,8 @@ subroutine DustEmissionK14( km, t_soil, w_top, rho_air,    &
    i2 = ubound(t_soil,1)
    j2 = ubound(t_soil,2)
 
-   allocate(w_g(i2,j2), w_gt(i2,j2), f_veg(i2,j2), clay(i2,j2), silt(i2,j2), k_gamma(i2,j2))
-   allocate(z0s(i2,j2), Dp_size(i2,j2))
+   allocate(w_g(i2,j2), w_gt(i2,j2), f_veg(i2,j2), clay(i2,j2), silt(i2,j2), k_gamma(i2,j2), source=0.0)
+   allocate(z0s(i2,j2), Dp_size(i2,j2), source=0.0)
 
    ! typical size of soil particles for optimal saltation is about 75e-6m
    Dp_size = 75e-6
@@ -1183,10 +1183,8 @@ subroutine Chem_SettlingSimple ( km, klid, flag, cdt, grav, &
 
    hsurf => hghte(i1:i2,j1:j2,km)
 
-   allocate(dz(i2,j2,km), radius(i2,j2,km), rhop(i2,j2,km), vsettle(i2,j2,km), qa(i2,j2,km))
-   allocate(cmass_before(i2,j2), cmass_after(i2,j2))
-   cmass_before = 0.d0
-   cmass_after = 0.d0
+   allocate(dz(i2,j2,km), radius(i2,j2,km), rhop(i2,j2,km), vsettle(i2,j2,km), qa(i2,j2,km), source=0.0)
+   allocate(cmass_before(i2,j2), cmass_after(i2,j2), source=0.0_DP)
 
    qa = int_qa
 
@@ -1359,10 +1357,8 @@ subroutine Chem_Settling ( km, klid, bin, flag, cdt, grav, &
 
    hsurf => hghte(i1:i2,j1:j2,km)
 
-   allocate(dz(i2,j2,km), radius(i2,j2,km), rhop(i2,j2,km), vsettle(i2,j2,km), qa(i2,j2,km))
-   allocate(cmass_before(i2,j2), cmass_after(i2,j2))
-   cmass_before = 0.d0
-   cmass_after = 0.d0
+   allocate(dz(i2,j2,km), radius(i2,j2,km), rhop(i2,j2,km), vsettle(i2,j2,km), qa(i2,j2,km), source=0.0)
+   allocate(cmass_before(i2,j2), cmass_after(i2,j2), source=0.0_DP)
 
    qa = int_qa
 
@@ -1967,7 +1963,7 @@ subroutine Chem_Settling2Gorig (km, klid, flag, bin, int_qa, grav, delp, &
    real, dimension(:,:,:), allocatable  :: vsettle   ! fall speed [m s-1]
    real(kind=DP), dimension(:,:,:), allocatable   :: dzd, vsd, qa, qa_temp
    real(kind=DP), dimension(:,:), allocatable  :: cmass_before, cmass_after
-   real(kind=DP) :: qdel, qsrc, d_p, dpm1
+   real(kind=DP) :: qdel, qsrc, d_p, d_pm1
 
    integer :: status
 
@@ -1986,9 +1982,10 @@ subroutine Chem_Settling2Gorig (km, klid, flag, bin, int_qa, grav, delp, &
 
 !  Allocate arrays
 !  ---------------
-   allocate(dz, mold=rhoa);
-   allocate(dzd(i2,j2,km), vsd(i2,j2,km), qa(i2,j2,km), vsettle(i2,j2,km), qa_temp(i2,j2,km))
-   allocate(cmass_before(i2,j2), cmass_after(i2,j2))
+   allocate(dz, mold=rhoa)
+   allocate(vsettle(i2,j2,km), source=0.0)
+   allocate(dzd(i2,j2,km), vsd(i2,j2,km), qa(i2,j2,km), qa_temp(i2,j2,km), source=0.0_DP)
+   allocate(cmass_before(i2,j2), cmass_after(i2,j2), source=0.0_DP)
 
 !  Handle the fact that hghte may be in the range [1,km+1] or [0,km]
 !  -----------------------------------------------------------------
@@ -2119,8 +2116,8 @@ subroutine Chem_Settling2Gorig (km, klid, flag, bin, int_qa, grav, delp, &
 !             do k = 2, km
              do k = klid+1, km
                d_p  = delp(i,j,k)
-               dpm1 = delp(i,j,k-1)
-               qsrc = qdel * dpm1 / d_p
+               d_pm1 = delp(i,j,k-1)
+               qsrc = qdel * d_pm1 / d_p
                qdel = qa(i,j,k)*dt_settle*vsd(i,j,k)/dzd(i,j,k)
                qa(i,j,k) = qa(i,j,k) - qdel + qsrc
             end do
@@ -2313,7 +2310,7 @@ subroutine Chem_SettlingSimpleOrig ( km, klid, flag, grav, cdt, radiusInp, rhopI
    real, dimension(:,:,:), allocatable  :: vsettle   ! fall speed [m s-1]
    real(kind=DP), dimension(:,:,:), allocatable   :: dzd, vsd, qa, qa_temp
    real(kind=DP), dimension(:,:), allocatable  :: cmass_before, cmass_after, qdel, &
-        dp, dpm1, qsrc
+        d_p, d_pm1, qsrc
 
 !EOP
 !-------------------------------------------------------------------------
@@ -2328,10 +2325,11 @@ subroutine Chem_SettlingSimpleOrig ( km, klid, flag, grav, cdt, radiusInp, rhopI
 
 !  Allocate arrays
 !  ---------------
-   allocate(dz, mold=rhoa);
-   allocate(dzd(i2,j2,km), vsd(i2,j2,km), qa(i2,j2,km), vsettle(i2,j2,km), qa_temp(i2,j2,km))
-   allocate(cmass_before(i2,j2), cmass_after(i2,j2), qdel(i2,j2), dp(i2,j2), &
-            dpm1(i2,j2), qsrc(i2,j2))
+   allocate(dz, mold=rhoa)
+   allocate(vsettle(i2,j2,km), source=0.0)
+   allocate(dzd(i2,j2,km), vsd(i2,j2,km), qa(i2,j2,km), qa_temp(i2,j2,km), source=0.0_DP)
+   allocate(cmass_before(i2,j2), cmass_after(i2,j2), qdel(i2,j2), d_p(i2,j2), &
+            d_pm1(i2,j2), qsrc(i2,j2), source=0.0_DP)
 
 !  Handle the fact that hghte may be in the range [1,km+1] or [0,km]
 !  -----------------------------------------------------------------
@@ -2459,9 +2457,9 @@ subroutine Chem_SettlingSimpleOrig ( km, klid, flag, grav, cdt, radiusInp, rhopI
      qa(i1:i2,j1:j2,klid) = qa(i1:i2,j1:j2,klid) - qdel
 
      do k = klid+1, km
-      dp   = delp(i1:i2,j1:j2,k)
-      dpm1 = delp(i1:i2,j1:j2,k-1)
-      qsrc = qdel * dpm1 / dp
+      d_p   = delp(i1:i2,j1:j2,k)
+      d_pm1 = delp(i1:i2,j1:j2,k-1)
+      qsrc = qdel * d_pm1 / d_p
       qdel = qa(i1:i2,j1:j2,k)*dt_settle*vsd(i1:i2,j1:j2,k)/dzd(i1:i2,j1:j2,k)
       qa(i1:i2,j1:j2,k) = qa(i1:i2,j1:j2,k) - qdel + qsrc
      enddo
@@ -2569,7 +2567,7 @@ subroutine DryDeposition ( km, tmpu, rhoa, hghte, oro, ustar, pblh, shflux, &
    i2 = dims(1); j2 = dims(2)
 
    allocate(dz(i2,j2),rmu(i2,j2),Ra(i2,j2),Rs(i2,j2),vdep(i2,j2), &
-            obk(i2,j2))
+            obk(i2,j2), source=0.0)
 
 !  Calculate the viscosity and thickness of the surface level
    dz = hghte(:,:,km-1) - hghte(:,:,km)
@@ -2818,7 +2816,7 @@ subroutine WetRemovalGOCART2G ( km, klid, n1, n2, bin_ind, cdt, aero_type, kin,
 
 !  Allocate arrays
    allocate(c_h2o(i2,j2,km), cldliq(i2,j2,km), cldice(i2,j2,km), pdog(i2,j2,km), &
-            delz(i2,j2,km), dpfli(i2,j2,km))
+            delz(i2,j2,km), dpfli(i2,j2,km), source=0.0)
 
 !  Initialize local variables
 !  --------------------------
@@ -2836,8 +2834,8 @@ subroutine WetRemovalGOCART2G ( km, klid, n1, n2, bin_ind, cdt, aero_type, kin,
    Td_ls = cdt
    nbins = n2-n1+1
 
-   allocate(fd(km,nbins),stat=ios)
-   allocate(dc(nbins),stat=ios)
+   allocate(fd(km,nbins),source=0.0,stat=ios)
+   allocate(dc(nbins),source=0.0,stat=ios)
 
    if( associated(fluxout) ) fluxout(i1:i2,j1:j2,bin_ind) = 0.0
 
@@ -3350,8 +3348,8 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch
    nch = size(channels)
    i2 = size(rhoa,1)
    j2 = size(rhoa,2)
-   allocate(fPMfm(nbins))
-   allocate(fPM25(nbins))
+   allocate(fPMfm(nbins),source=0.0)
+   allocate(fPM25(nbins),source=0.0)
 
 !  Get the wavelength indices
 !  --------------------------
@@ -3416,7 +3414,7 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch
       ilam470 .ne. ilam870) do_angstrom = .true.
 
    if( present(angstrom) .and. do_angstrom ) then
-      allocate(tau470(i1:i2,j1:j2), tau870(i1:i2,j1:j2))
+      allocate(tau470(i1:i2,j1:j2), tau870(i1:i2,j1:j2), source=0.0)
    end if
 
 !  Compute the fine mode (sub-micron) and PM2.5 bin-wise fractions
@@ -4240,6 +4238,7 @@ subroutine hoppelCorrection (radius, rhop, rh, dz, ustar, rhFlag, &
    rc = __SUCCESS__
    fhoppel = 1.0
    allocate(vsettle, mold=rh)
+   vsettle=0.0
 
    do j = 1, ubound(rh,2)
       do i = 1, ubound(rh,1)
@@ -4376,7 +4375,7 @@ subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraf
    ijl = ( i2 - i1 + 1 ) * ( j2 - j1 + 1 )
 
    allocate(factor(i2,j2), srcHydrophobic(i2,j2), srcHydrophilic(i2,j2), srcBiofuel(i2,j2), &
-            srcBiomass(i2,j2), srcAnthro(i2,j2), srcBiogenic(i2,j2), f_bb_(i2,j2), exttau_bb_(i2,j2))
+            srcBiomass(i2,j2), srcAnthro(i2,j2), srcBiogenic(i2,j2), f_bb_(i2,j2), exttau_bb_(i2,j2), source=0.0)
 
 !  Emission factors scaling from source files to desired mass quantity
    eBiomass = ratPOM
@@ -4427,6 +4426,11 @@ subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraf
    allocate(p500, mold=pblh)
    allocate(pPBL, mold=pblh)
    ps = 0.0
+   p0 = 0.0
+   z0 = 0.0
+   p100 = 0.0
+   p500 = 0.0
+   pPBL = 0.0
    do k = 1, km
     ps(i1:i2,j1:j2) = ps(i1:i2,j1:j2) + delp(i1:i2,j1:j2,k)
    end do
@@ -5454,7 +5458,7 @@ subroutine SulfateDistributeEmissions ( km, nbins, cdt, grav, nymd, nhms, &
    hsurf = hghte(i1:i2,j1:j2,km)
 
    allocate(srcSO2(i2,j2), srcSO4(i2,j2), srcDMS(i2,j2), srcSO4anthro(i2,j2), &
-            srcSO2anthro(i2,j2), srcSO2bioburn(i2,j2))
+            srcSO2anthro(i2,j2), srcSO2bioburn(i2,j2), source=0.0)
 
 !  Initialize local variables
 !  --------------------------
@@ -5505,6 +5509,11 @@ subroutine SulfateDistributeEmissions ( km, nbins, cdt, grav, nymd, nhms, &
    allocate(pPblh, mold=pblh)
 
    ps = 0.0
+   p0 = 0.0
+   z0 = 0.0
+   p100 = 0.0
+   p500 = 0.0
+   pPblh = 0.0
    do k = 1, km
     ps(i1:i2,j1:j2) = ps(i1:i2,j1:j2) + delp(i1:i2,j1:j2,k)
    end do
@@ -5942,7 +5951,7 @@ subroutine SulfateUpdateOxidants (nymd_current, nhms_current, lonRad, latRad, &
     j2 = size(rhoa,2)
 
     allocate(cossza(i1:i2,j1:j2), sza(i1:i2,j1:j2), tcosz(i1:i2,j1:j2), &
-             tday(i1:i2,j1:j2), tnight(i1:i2,j1:j2))
+             tday(i1:i2,j1:j2), tnight(i1:i2,j1:j2), source=0.0)
 
 ! Update emissions/production if necessary (daily)
 !  -----------------------------------------------
@@ -6845,7 +6854,7 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t
    j2 = ubound(tmpu, 2)
    i2 = ubound(tmpu, 1)
 
-   allocate(tau470(i1:i2,j1:j2), tau870(i1:i2,j1:j2))
+   allocate(tau470(i1:i2,j1:j2), tau870(i1:i2,j1:j2), source=0.0)
 
 !  Get the wavelength indices
 !  --------------------------
@@ -7263,6 +7272,10 @@ subroutine SulfateChemDriver (km, klid, cdt, PI, radToDeg, von_karman, &
    allocate(cossza, mold=oro)
    allocate(sza, mold=oro)
 
+   drydepositionfrequency = 0.0
+   cossza = 0.0
+   sza = 0.0
+
 !  Reset the production terms
    allocate(pSO2_DMS, mold=tmpu)
    allocate(pMSA_DMS, mold=tmpu)
@@ -7464,6 +7477,8 @@ subroutine SulfateChemDriver_DMS (km, klid, cdt, airMolWght, nAvogadro, cpd, &
 
    allocate(pSO2_DMS, mold=tmpu)
    allocate(pMSA_DMS, mold=tmpu)
+   pSO2_DMS = 0.0
+   pMSA_DMS = 0.0
 
 !  spatial loop
    do k = klid, km
@@ -7629,6 +7644,9 @@ subroutine SulfateChemDriver_SO2 (km, klid, cdt, airMolWght, nAvogadro, cpd, gra
    allocate(pSO4g_SO2, mold=tmpu)
    allocate(pSO4aq_SO2, mold=tmpu)
    allocate(fout(i2,j2))
+   pSO4g_SO2 = 0.0
+   pSO4aq_SO2 = 0.0
+   fout = 0.0
 
 !  Conversion of SO2 mmr to SO2 vmr
    fMR = airMolWght / fMassSO2
@@ -7895,6 +7913,7 @@ subroutine SulfateChemDriver_MSA (km, klid, cdt, grav, qa, nMSA, delp, drydepf,
    i2 = ubound(qa, 1)
 
    allocate(fout(i2,j2))
+   fout = 0.0
 
 !  spatial loop
    do k = klid, km