Implement bromoform tracer and coupling

cesm/mod_cesm.F90

@@ -28,7 +28,8 @@ module mod_cesm
    use mod_time, only: nstep
    use mod_xc
    use mod_forcing, only: trxday, srxday, swa, nsf, lip, sop, eva, rnf, rfi, &
-                          fmltfz, sfl, ztx, mty, ustarw, slp, abswnd, atmco2
+                          fmltfz, sfl, ztx, mty, ustarw, slp, abswnd, atmco2,&
+                          atmbrf
    use mod_ben02, only: initai, rdcsic, rdctsf, fnlzai
    use mod_seaice, only: ficem
    use mod_checksum, only: csdiag, chksummsk
@@ -70,6 +71,7 @@ module mod_cesm
       slp_da, &          ! Sea-level pressure [kg m-1 s-2].
       abswnd_da, &       ! Wind speed at measurement height (zu) [m s-1].
       atmco2_da, &       ! Atmospheric CO2 concentration [ppm].
+      atmbrf_da, &       ! Atmospheric bromoform concentration [ppt].
       ficem_da           ! Ice concentration [].
 
    logical :: &
@@ -81,7 +83,7 @@ module mod_cesm
    public :: runid_cesm, runtyp_cesm, ocn_cpl_dt_cesm, nstep_in_cpl, hmlt, &
              frzpot, mltpot, swa_da, nsf_da, hmlt_da, lip_da, sop_da, eva_da, &
              rnf_da, rfi_da, fmltfz_da, sfl_da, ztx_da, mty_da, ustarw_da, &
-             slp_da, abswnd_da, atmco2_da, ficem_da, smtfrc, l1ci, l2ci, &
+             slp_da, abswnd_da, atmco2_da, atmbrf_da, ficem_da, smtfrc, l1ci, l2ci, &
              inicon_cesm, inifrc_cesm, getfrc_cesm
 
 contains
@@ -176,6 +178,7 @@ subroutine getfrc_cesm
            ficem(i, j)  = w1*ficem_da(i, j, l1ci)  + w2*ficem_da(i, j, l2ci)
            abswnd(i, j) = w1*abswnd_da(i, j, l1ci) + w2*abswnd_da(i, j, l2ci)
            atmco2(i, j) = w1*atmco2_da(i, j, l1ci) + w2*atmco2_da(i, j, l2ci)
+           atmbrf(i, j) = w1*atmbrf_da(i, j, l1ci) + w2*atmbrf_da(i, j, l2ci)
         enddo
         enddo
         do l = 1, isu(j)
@@ -210,6 +213,7 @@ subroutine getfrc_cesm
       call ncdefvar('ficem_da', 'x y', ndouble, 8)
       call ncdefvar('abswnd_da', 'x y', ndouble, 8)
       call ncdefvar('atmco2_da', 'x y', ndouble, 8)
+      call ncdefvar('atmbrf_da', 'x y', ndouble, 8)
       call ncdefvar('ztx_da', 'x y', ndouble, 8)
       call ncdefvar('mty_da', 'x y', ndouble, 8)
       call ncedef
@@ -244,6 +248,8 @@ subroutine getfrc_cesm
                   ip, 1, 1._r8, 0._r8, 8)
       call ncwrtr('atmco2_da', 'x y', atmco2_da(1 - nbdy, 1 - nbdy, l2ci), &
                   ip, 1, 1._r8, 0._r8, 8)
+      call ncwrtr('atmbrf_da', 'x y', atmbrf_da(1 - nbdy, 1 - nbdy, l2ci), &
+     .            ip, 1, 1._r8, 0._r8, 8)
       call ncwrtr('ztx_da', 'x y', ztx_da(1 - nbdy, 1 - nbdy, l2ci), &
                   iu, 1, 1._r8, 0._r8, 8)
       call ncwrtr('mty_da', 'x y', mty_da(1 - nbdy, 1 - nbdy, l2ci), &
@@ -274,6 +280,7 @@ subroutine getfrc_cesm
          call chksummsk(ficem, ip, 1, 'ficem')
          call chksummsk(abswnd, ip, 1, 'abswnd')
          call chksummsk(atmco2, ip, 1, 'atmco2')
+         call chksummsk(atmbrf, ip, 1, 'atmbrf')
       endif
 
    end subroutine getfrc_cesm

cime_config/buildlib

@@ -21,6 +21,7 @@ set HAMOCC_CFC = `./xmlquery HAMOCC_CFC --value`
 set HAMOCC_NATTRC = `./xmlquery HAMOCC_NATTRC --value`
 set HAMOCC_SEDBYPASS = `./xmlquery HAMOCC_SEDBYPASS --value`
 set HAMOCC_CISO = `./xmlquery HAMOCC_CISO --value`
+set HAMOCC_VSLS = `./xmlquery HAMOCC_VSLS --value`
 set GMAKE = `./xmlquery GMAKE --value`
 set GMAKE_J = `./xmlquery  GMAKE_J --value`
 set OCN_CO2_TYPE = `./xmlquery OCN_CO2_TYPE --value`
@@ -109,6 +110,9 @@ if ($#tracers != 0) then
         endif
         set cpp_ocn = "$cpp_ocn -Dcisonew"
       endif       
+      if ($HAMOCC_VSLS == TRUE) then
+        set cpp_ocn = "$cpp_ocn -DBROMO"
+      endif       
       if ($co2type == prognostic) then
         set cpp_ocn = "$cpp_ocn -DPROGCO2"
       else if ($co2type == diagnostic) then

cime_config/buildnml

@@ -408,6 +408,8 @@ set SRF_CO214FXU  = '0, 2, 2'
 set SRF_CFC11     = '0, 2, 2' 
 set SRF_CFC12     = '0, 2, 2' 
 set SRF_SF6       = '0, 2, 2' 
+set SRF_BROMO     = '0, 2, 2'
+set SRF_BROMOFX   = '0, 2, 2'
 set INT_PHOSY     = '4, 2, 2' 
 set INT_NFIX      = '0, 2, 2' 
 set INT_DNIT      = '0, 2, 2' 
@@ -479,6 +481,7 @@ set LYR_WPHY      = '0, 0, 2'
 set LYR_WNOS      = '0, 0, 2'
 set LYR_EPS       = '0, 0, 0'
 set LYR_ASIZE     = '0, 0, 0'
+set LYR_BROMO     = '0, 0, 2'
 set BGC_DP        = '0, 2, 2'
 set LVL_PHYTO     = '0, 2, 2'
 set LVL_GRAZER    = '0, 2, 2'
@@ -528,6 +531,7 @@ set LVL_WPHY      = '0, 2, 2'
 set LVL_WNOS      = '0, 2, 2'
 set LVL_EPS       = '0, 0, 0'
 set LVL_ASIZE     = '0, 0, 0'
+set LVL_BROMO     = '0, 2, 2'
 set SDM_POWAIC    = '0, 0, 2'
 set SDM_POWAAL    = '0, 0, 2'
 set SDM_POWAPH    = '0, 0, 2'
@@ -1502,6 +1506,8 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
   SRF_CFC11     = $SRF_CFC11
   SRF_CFC12     = $SRF_CFC12
   SRF_SF6       = $SRF_SF6
+  SRF_BROMO     = $SRF_BROMO
+  SRF_BROMOFX   = $SRF_BROMOFX
   INT_PHOSY     = $INT_PHOSY
   INT_NFIX      = $INT_NFIX
   INT_DNIT      = $INT_DNIT
@@ -1574,6 +1580,7 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
   LYR_EPS       = $LYR_EPS
   LYR_ASIZE     = $LYR_ASIZE
   LYR_DP        = $BGC_DP
+  LYR_BROMO     = $LYR_BROMO
   LVL_PHYTO     = $LVL_PHYTO
   LVL_GRAZER    = $LVL_GRAZER
   LVL_DOC       = $LVL_DOC
@@ -1622,6 +1629,7 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
   LVL_WNOS      = $LVL_WNOS
   LVL_EPS       = $LVL_EPS
   LVL_ASIZE     = $LVL_ASIZE
+  LVL_BROMO     = $LVL_BROMO
   SDM_POWAIC    = $SDM_POWAIC
   SDM_POWAAL    = $SDM_POWAAL
   SDM_POWAPH    = $SDM_POWAPH

cime_config/config_component.xml

@@ -147,6 +147,18 @@
     <file>env_build.xml</file>
     <desc>Set preprocessor option to activate the carbon isotope code. Requires module ecosys</desc>
   </entry>
+
+  <entry id="HAMOCC_VSLS">
+    <type>logical</type>
+    <valid_values>TRUE,FALSE</valid_values>
+    <default_value>FALSE</default_value>
+    <values>
+      <value compset="%VSLSC">TRUE</value>
+    </values>
+    <group>build_component_blom</group>
+    <file>env_build.xml</file>
+    <desc>Set preprocessor option to activate the VSLS-Bromoform tracer code. Requires module ecosys</desc>
+  </entry>
 
   <entry id="BLOM_TURBULENT_CLOSURE">
     <type>char</type>

drivers/cpl_mct/export_mct.F

@@ -148,6 +148,24 @@ subroutine export_mct(o2x_o, lsize, perm, jjcpl, nsend, sbuff,
      .      write (lp,*) 'export_mct: co2 flux not sent to coupler'
       endif
 
+      ! ----------------------------------------------------------------
+      ! Pack bromoform flux (kg CHBr3/m^2/s), if requested
+      ! ----------------------------------------------------------------
+
+      if (index_o2x_Faoo_fbrf_ocn > 0) then
+        n = 0
+         do j = 1, jjcpl
+            do i = 1, ii
+               n = n + 1
+               o2x_o%rAttr(index_o2x_Faoo_fbrf_ocn,n) =
+     .            sbuff(i,j,index_o2x_Faoo_fbrf_ocn)*tfac
+            enddo
+         enddo
+      else
+         if (mnproc.eq.1)
+     .     write (lp,*) 'export_mct: bromoform flux not sent to coupler'
+      endif
+
       tlast_coupled = 0._r8
 
       !-----------------------------------------------------------------

drivers/cpl_mct/import_mct.F

@@ -29,7 +29,7 @@ subroutine import_mct(x2o_o, lsize, perm, jjcpl)
       use mod_cesm, only: swa_da, nsf_da, hmlt_da, lip_da, sop_da,
      .                    eva_da, rnf_da, rfi_da, fmltfz_da, sfl_da,
      .                    ztx_da, mty_da, ustarw_da, slp_da, abswnd_da,
-     .                    atmco2_da, ficem_da, l1ci, l2ci
+     .                    atmco2_da, atmbrf_da, ficem_da, l1ci, l2ci
       use mod_utility, only: util1, util2
       use mod_checksum, only: csdiag, chksummsk
       use blom_cpl_indices
@@ -309,6 +309,40 @@ subroutine import_mct(x2o_o, lsize, perm, jjcpl)
      .   write (lp,*) 'import_mct: atmospheric co2 not read'
 #endif      
 
+      if (index_x2o_Sa_brfprog > 0) then
+         n = 0
+         do j = 1, jjcpl
+            do i = 1, ii
+               n = n + 1
+               if     (ip(i,j) == 0) then
+                  atmbrf_da(i,j,l2ci) = mval
+               elseif (cplmsk(i,j) == 0) then
+                  atmbrf_da(i,j,l2ci) = fval
+               else
+                  ! Atmospheric bromoform concentration [ppt]
+                  atmbrf_da(i,j,l2ci) =
+     .               x2o_o%rAttr(index_x2o_Sa_brfprog,n)
+               endif
+            enddo
+         enddo
+         call fill_global(mval, fval, halo_ps,
+     .                    atmbrf_da(1-nbdy,1-nbdy,l2ci))
+         if (mnproc.eq.1)
+     .      write (lp,*) 'import_mct: prog. atmospheric bromoform read'
+      else
+         do j = 1, jj
+            do i = 1, ii
+               if (ip(i,j) == 0) then
+                  atmbrf_da(i,j,l2ci) = mval
+               else
+                  atmbrf_da(i,j,l2ci) = -1
+               endif
+            enddo
+         enddo
+         if (mnproc.eq.1)
+     .   write (lp,*) 'import_mct: prog. atmospheric bromoform not read'
+      endif
+
       if (csdiag) then
          if (mnproc.eq.1) then
             write (lp,*) 'import_mct:'
@@ -329,7 +363,8 @@ subroutine import_mct(x2o_o, lsize, perm, jjcpl)
          call chksummsk(slp_da(1-nbdy,1-nbdy,l2ci),ip,1,'slp')
          call chksummsk(ficem_da(1-nbdy,1-nbdy,l2ci),ip,1,'ficem')
          call chksummsk(abswnd_da(1-nbdy,1-nbdy,l2ci),ip,1,'abswnd')
-         call chksummsk(atmco2_da(1-nbdy,1-nbdy,l2ci),ip,1,'abswnd')
+         call chksummsk(atmco2_da(1-nbdy,1-nbdy,l2ci),ip,1,'atmco2')
+         call chksummsk(atmbrf_da(1-nbdy,1-nbdy,l2ci),ip,1,'atmbrf')
       endif
 
       end subroutine import_mct

drivers/cpl_mct/sumsbuff_mct.F

@@ -27,7 +27,7 @@ subroutine sumsbuff_mct(nsend, sbuff, tlast_coupled)
       use mod_grid, only: scuy, scvx, scuxi, scvyi
       use mod_state, only: u, v, temp, saln, pbu, pbv, ubflxs, vbflxs,
      .                     sealv
-      use mod_forcing, only: flxco2, flxdms
+      use mod_forcing, only: flxco2, flxdms, flxbrf
       use mod_cesm, only: frzpot
       use blom_cpl_indices
 
@@ -141,6 +141,18 @@ subroutine sumsbuff_mct(nsend, sbuff, tlast_coupled)
          enddo
       endif
 
+      if (index_o2x_Faoo_fbrf_ocn > 0) then
+         do j = 1, jj
+            do l = 1, isp(j)
+            do i = max(1,ifp(j,l)), min(ii,ilp(j,l))
+               sbuff(i,j,index_o2x_Faoo_fbrf_ocn) =
+     .            sbuff(i,j,index_o2x_Faoo_fbrf_ocn)
+     .          + flxbrf(i,j)*baclin
+            enddo
+            enddo
+         enddo
+      endif
+
 
       !-----------------------------------------------------------------
       ! Increment time since last coupling

drivers/cpl_share/blom_cpl_indices.F90

@@ -43,6 +43,7 @@ module blom_cpl_indices
   integer :: index_o2x_Fioo_q
   integer :: index_o2x_Faoo_fco2_ocn
   integer :: index_o2x_Faoo_fdms_ocn
+  integer :: index_o2x_Faoo_fbrf_ocn
 
   ! drv -> ocn
 
@@ -51,6 +52,7 @@ module blom_cpl_indices
   integer :: index_x2o_Sa_pslv         ! sea-level pressure               (Pa)
   integer :: index_x2o_Sa_co2prog      ! bottom atm level prognostic CO2
   integer :: index_x2o_Sa_co2diag      ! bottom atm level diagnostic CO2
+  integer :: index_x2o_Sa_brfprog      ! bottom atm level prognostic bromoform (ppt)
   integer :: index_x2o_Faxa_nhx        ! nitrogen deposition (nhx) flux from atm (kgNm2/sec)
   integer :: index_x2o_Faxa_noy        ! nitrogen deposition (noy) flux from atm (kgNm2/sec)
 
@@ -126,6 +128,7 @@ subroutine blom_cpl_indices_set( )
     index_o2x_Fioo_q        = mct_avect_indexra(o2x,'Fioo_q')
     index_o2x_Faoo_fco2_ocn = mct_avect_indexra(o2x,'Faoo_fco2_ocn',perrWith='quiet')
     index_o2x_Faoo_fdms_ocn = mct_avect_indexra(o2x,'Faoo_fdms_ocn',perrWith='quiet')
+    index_o2x_Faoo_fbrf_ocn = mct_avect_indexra(o2x,'Faoo_fbrf_ocn',perrWith='quiet')
     index_x2o_Si_ifrac      = mct_avect_indexra(x2o,'Si_ifrac')
     index_x2o_Sa_pslv       = mct_avect_indexra(x2o,'Sa_pslv')
     index_x2o_So_duu10n     = mct_avect_indexra(x2o,'So_duu10n')
@@ -169,6 +172,7 @@ subroutine blom_cpl_indices_set( )
     index_x2o_Faxa_dstwet4  = mct_avect_indexra(x2o,'Faxa_dstwet4')
     index_x2o_Sa_co2prog    = mct_avect_indexra(x2o,'Sa_co2prog',perrWith='quiet')
     index_x2o_Sa_co2diag    = mct_avect_indexra(x2o,'Sa_co2diag',perrWith='quiet')
+    index_x2o_Sa_brfprog    = mct_avect_indexra(x2o,'Sa_brfprog',perrWith='quiet')
     index_x2o_Faxa_nhx      = mct_avect_indexra(x2o,'Faxa_nhx',perrWith='quiet')
     index_x2o_Faxa_noy      = mct_avect_indexra(x2o,'Faxa_noy',perrWith='quiet')
 

hamocc/accfields.F90

@@ -127,6 +127,10 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
 #ifdef natDIC
       call accsrf(jnatco2fx,atmflx(1,1,iatmnco2),omask,0)
 #endif
+#ifdef BROMO
+      call accsrf(jatmbromo,atm(1,1,iatmbromo),omask,0)
+      call accsrf(jbromofx,atmflx(1,1,iatmbromo),omask,0)
+#endif
 #ifdef cisonew
       call accsrf(jatmc13,atm(1,1,iatmc13),omask,0)
       call accsrf(jatmc14,atm(1,1,iatmc14),omask,0)
@@ -168,6 +172,9 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
       call accsrf(jsrfnatalk,ocetra(1,1,1,inatalkali),omask,0)
       call accsrf(jnatpco2,natpco2d,omask,0)
 #endif
+#ifdef BROMO
+      call accsrf(jsrfbromo,ocetra(1,1,1,ibromo),omask,0)
+#endif
 
 ! Accumulate the diagnostic mass sinking field 
       IF( domassfluxes ) THEN
@@ -251,6 +258,9 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
       call acclyr(jcfc12,ocetra(1,1,1,icfc12),pddpo,1)
       call acclyr(jsf6,ocetra(1,1,1,isf6),pddpo,1)
 #endif
+#ifdef BROMO
+      call acclyr(jbromo,ocetra(1,1,1,ibromo),pddpo,1)
+#endif
 
 
 ! Accumulate level diagnostics
@@ -262,7 +272,7 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
      &  jlvlnatomegaa+jlvlnatomegac+jlvldic13+jlvldic14+jlvld13c+       &
      &  jlvld14c+jlvlbigd14c+jlvlpoc13+jlvldoc13+jlvlcalc13+jlvlphyto13+&
      &  jlvlgrazer13+jlvlnos+jlvlwphy+jlvlwnos+jlvleps+jlvlasize+       &
-     &  jlvlcfc11+jlvlcfc12+jlvlsf6).NE.0) THEN
+     &  jlvlcfc11+jlvlcfc12+jlvlsf6+jlvlbromo).NE.0) THEN
         DO k=1,kpke
           call bgczlv(pddpo,k,ind1,ind2,wghts)
           call acclvl(jlvlphyto,ocetra(1,1,1,iphy),k,ind1,ind2,wghts)
@@ -322,6 +332,9 @@ SUBROUTINE ACCFIELDS(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask)
           call acclvl(jlvlcfc11,ocetra(1,1,1,icfc11),k,ind1,ind2,wghts)
           call acclvl(jlvlcfc12,ocetra(1,1,1,icfc12),k,ind1,ind2,wghts)
           call acclvl(jlvlsf6,ocetra(1,1,1,isf6),k,ind1,ind2,wghts)
+#endif
+#ifdef BROMO
+          call acclvl(jlvlbromo,ocetra(1,1,1,ibromo),k,ind1,ind2,wghts)
 #endif
         ENDDO
       ENDIF

hamocc/aufr_bgc.F90

@@ -129,7 +129,7 @@ SUBROUTINE AUFR_BGC(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,               &
       INTEGER   :: restdtoce                         !  time step number from bgc ocean file
       INTEGER   :: idate(5),i,j,k
       character :: rstfnm*256
-      logical   :: lread_cfc,lread_nat,lread_iso,lread_atm
+      logical   :: lread_cfc,lread_nat,lread_iso,lread_atm,lread_bro
 #ifdef cisonew
       REAL :: rco213,rco214,alpha14,beta13,beta14,d13C_atm,d14cat
 #endif
@@ -349,6 +349,26 @@ SUBROUTINE AUFR_BGC(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,               &
       ENDIF
 #endif
 
+! Find out whether to restart Bromoform
+#ifdef BROMO
+      lread_bro=.true.
+      IF(IOTYPE==0) THEN
+        if(mnproc==1) ncstat=nf90_inq_varid(ncid,'bromo',ncvarid)
+        call xcbcst(ncstat)
+        if(ncstat.ne.nf90_noerr) lread_bro=.false.
+      ELSE IF(IOTYPE==1) THEN
+#ifdef PNETCDF
+        ncstat=nfmpi_inq_varid(ncid,'bromo',ncvarid)
+        if(ncstat.ne.nf_noerr) lread_bro=.false.
+#endif
+      ENDIF
+      IF(mnproc==1 .and. .not. lread_bro) THEN
+        WRITE(io_stdo_bgc,*) ' '
+        WRITE(io_stdo_bgc,*) 'AUFR_BGC info: Bromoform tracer not in restart file, '
+        WRITE(io_stdo_bgc,*) 'Initialised to 0.01 pmol L-1 (Stemmler et al., 2015).'
+      ENDIF
+#endif
+
 ! Find out whether to restart atmosphere
 #if defined(BOXATM)
       lread_atm=.true.
@@ -435,6 +455,12 @@ SUBROUTINE AUFR_BGC(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,               &
       CALL read_netcdf_var(ncid,'hi',nathi(1,1,1),kpke,0,iotype)
       ENDIF
 #endif
+#ifdef BROMO
+      IF(lread_bro) THEN
+      CALL read_netcdf_var(ncid,'bromo',locetra(1,1,1,ibromo),2*kpke,0,iotype)
+      ENDIF
+#endif
+
 !
 ! Read restart data : diagnostic ocean fields (needed for bit to bit reproducability)
 !