Add capability to run C-isotopes with sediment code switched on (test…

…ed technically but not validated, build script will still refuse to run this combination of options)
NorESMhub · Aug 3, 2023 · aa737aa · aa737aa
1 parent 4930ba4
commit aa737aa
Show file tree

Hide file tree

Showing 4 changed files with 79 additions and 23 deletions.
diff --git a/hamocc/aufr_bgc.F90 b/hamocc/aufr_bgc.F90
@@ -482,11 +482,14 @@ SUBROUTINE AUFR_BGC(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,               &
       CALL read_netcdf_var(ncid,'powasi',powtra2(1,1,1,ipowasi),2*ks,0,iotype)
 #ifdef cisonew
       IF(lread_iso) THEN
-      ! Burial fields for c-isotopes still missing
       CALL read_netcdf_var(ncid,'ssso13',sedlay2(1,1,1,issso13),2*ks,0,iotype)
       CALL read_netcdf_var(ncid,'ssso14',sedlay2(1,1,1,issso14),2*ks,0,iotype)
       CALL read_netcdf_var(ncid,'sssc13',sedlay2(1,1,1,isssc13),2*ks,0,iotype)
       CALL read_netcdf_var(ncid,'sssc14',sedlay2(1,1,1,isssc14),2*ks,0,iotype)
+      CALL read_netcdf_var(ncid,'bur_o13',burial2(1,1,1,issso13),2,0,iotype)
+      CALL read_netcdf_var(ncid,'bur_o14',burial2(1,1,1,issso14),2,0,iotype)
+      CALL read_netcdf_var(ncid,'bur_c13',burial2(1,1,1,isssc13),2,0,iotype)
+      CALL read_netcdf_var(ncid,'bur_c14',burial2(1,1,1,isssc14),2,0,iotype)
       CALL read_netcdf_var(ncid,'powc13',powtra2(1,1,1,ipowc13),2*ks,0,iotype)
       CALL read_netcdf_var(ncid,'powc14',powtra2(1,1,1,ipowc14),2*ks,0,iotype)
       ENDIF
@@ -586,18 +589,34 @@ SUBROUTINE AUFR_BGC(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,               &
       ENDDO
 #ifndef sedbypass
       ! Burial fields for c-isotopes still missing
+      !JT added burial loop below 20.06.2023
       DO  k=1,2*ks
       DO  j=1,kpje
       DO  i=1,kpie 
         IF(omask(i,j) .GT. 0.5) THEN
         rco213=ocetra(i,j,kbo(i,j),isco213)/(ocetra(i,j,kbo(i,j),isco212)+safediv)
         rco214=ocetra(i,j,kbo(i,j),isco214)/(ocetra(i,j,kbo(i,j),isco212)+safediv)
-        powtra2(i,j,k,ipowc13)=powtra2(i,j,k,ipowc)*rco213*bifr13
-        powtra2(i,j,k,ipowc14)=powtra2(i,j,k,ipowc)*rco214*bifr14
-        sedlay2(i,j,k,issso13)=sedlay2(i,j,k,issso)*rco213*bifr13
-        sedlay2(i,j,k,issso14)=sedlay2(i,j,k,issso)*rco214*bifr14
-        sedlay2(i,j,k,isssc13)=sedlay2(i,j,k,isssc)*rco213
-        sedlay2(i,j,k,isssc14)=sedlay2(i,j,k,isssc)*rco214
+        powtra2(i,j,k,ipowc13)=powtra2(i,j,k,ipowaic)*rco213
+        powtra2(i,j,k,ipowc14)=powtra2(i,j,k,ipowaic)*rco214
+        sedlay2(i,j,k,issso13)=sedlay2(i,j,k,issso12)*rco213*bifr13
+        sedlay2(i,j,k,issso14)=sedlay2(i,j,k,issso12)*rco214*bifr14
+        sedlay2(i,j,k,isssc13)=sedlay2(i,j,k,isssc12)*rco213
+        sedlay2(i,j,k,isssc14)=sedlay2(i,j,k,isssc12)*rco214
+        ENDIF
+      ENDDO
+      ENDDO
+      ENDDO
+
+      DO  k=1,2
+      DO  j=1,kpje
+      DO  i=1,kpie 
+        IF(omask(i,j) .GT. 0.5) THEN
+        rco213=ocetra(i,j,kbo(i,j),isco213)/(ocetra(i,j,kbo(i,j),isco212)+safediv)
+        rco214=ocetra(i,j,kbo(i,j),isco214)/(ocetra(i,j,kbo(i,j),isco212)+safediv)
+        burial2(i,j,k,issso13)=burial2(i,j,k,issso12)*rco213*bifr13
+        burial2(i,j,k,issso14)=burial2(i,j,k,issso12)*rco214*bifr14
+        burial2(i,j,k,isssc13)=burial2(i,j,k,isssc12)*rco213
+        burial2(i,j,k,isssc14)=burial2(i,j,k,isssc12)*rco214
         ENDIF
       ENDDO
       ENDDO

diff --git a/hamocc/aufw_bgc.F90 b/hamocc/aufw_bgc.F90
@@ -733,6 +733,23 @@ SUBROUTINE AUFW_BGC(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,              &
      &    7,'kg/m**2',20,'Burial layer of clay',                        &
      &    rmissing,93,io_stdo_bgc)
 
+#ifdef cisonew
+      CALL NETCDF_DEF_VARDB(ncid,8,'bur_o13',3,ncdimst,ncvarid,         &
+     &    9,'kmol/m**2',27,'Burial layer of organic 13C',               &
+     &    rmissing,94,io_stdo_bgc)
+
+      CALL NETCDF_DEF_VARDB(ncid,8,'bur_o14',3,ncdimst,ncvarid,         &
+     &    9,'kmol/m**2',27,'Burial layer of organic 14C',               &
+     &    rmissing,95,io_stdo_bgc)
+
+      CALL NETCDF_DEF_VARDB(ncid,8,'bur_c13',3,ncdimst,ncvarid,         &
+     &    9,'kmol/m**2',23,'Burial layer of Ca13CO3',                   &
+     &    rmissing,96,io_stdo_bgc)
+
+      CALL NETCDF_DEF_VARDB(ncid,8,'bur_c14',3,ncdimst,ncvarid,         &
+     &    9,'kmol/m**2',23,'Burial layer of Ca14CO3',                   &
+     &    rmissing,97,io_stdo_bgc)
+#endif
 #endif /* sedbypass */
 
 !
@@ -897,6 +914,10 @@ SUBROUTINE AUFW_BGC(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,              &
       CALL write_netcdf_var(ncid,'ssso14',sedlay2(1,1,1,issso14),2*ks,0)
       CALL write_netcdf_var(ncid,'sssc13',sedlay2(1,1,1,isssc13),2*ks,0)
       CALL write_netcdf_var(ncid,'sssc14',sedlay2(1,1,1,isssc14),2*ks,0)
+      CALL write_netcdf_var(ncid,'bur_o13',burial2(1,1,1,issso13),2,0)
+      CALL write_netcdf_var(ncid,'bur_o14',burial2(1,1,1,issso14),2,0)
+      CALL write_netcdf_var(ncid,'bur_c13',burial2(1,1,1,isssc13),2,0)
+      CALL write_netcdf_var(ncid,'bur_c14',burial2(1,1,1,isssc14),2,0)
       CALL write_netcdf_var(ncid,'powc13',powtra2(1,1,1,ipowc13),2*ks,0)
       CALL write_netcdf_var(ncid,'powc14',powtra2(1,1,1,ipowc14),2*ks,0)
 #endif

diff --git a/hamocc/carchm.F90 b/hamocc/carchm.F90
@@ -637,6 +637,17 @@ SUBROUTINE CARCHM(kpie,kpje,kpke,kbnd,                                  &
         enddo
 !$OMP END PARALLEL DO
         enddo
+
+!$OMP PARALLEL DO PRIVATE(i)
+        do j=1,kpje
+        do i=1,kpie
+        if(omask(i,j).gt.0.5) then
+        burial(i,j,issso14) = burial(i,j,issso14)*c14dec
+        burial(i,j,isssc14) = burial(i,j,isssc14)*c14dec
+        endif
+        enddo
+        enddo
+!$OMP END PARALLEL DO
 #endif
 #endif
 

diff --git a/hamocc/sedshi.F90 b/hamocc/sedshi.F90
@@ -230,6 +230,16 @@ SUBROUTINE SEDSHI(kpie,kpje,omask)
         sedlay(i,j,ks,issster)=sedlay(i,j,ks,issster)                  &
      &                        +refill*burial(i,j,issster)/frac
 
+#ifdef cisonew
+        sedlay(i,j,ks,issso13)=sedlay(i,j,ks,issso13)                  &
+     &                        +refill*burial(i,j,issso13)/frac
+        sedlay(i,j,ks,isssc13)=sedlay(i,j,ks,isssc13)                  &
+     &                        +refill*burial(i,j,isssc13)/frac
+        sedlay(i,j,ks,issso14)=sedlay(i,j,ks,issso14)                  &
+     &                        +refill*burial(i,j,issso14)/frac
+        sedlay(i,j,ks,isssc14)=sedlay(i,j,ks,isssc14)                  &
+     &                        +refill*burial(i,j,isssc14)/frac
+#endif
 ! account for losses in buried sediment
         burial(i,j,issso12) = burial(i,j,issso12)                      &
      &                      - refill*burial(i,j,issso12)
@@ -239,7 +249,16 @@ SUBROUTINE SEDSHI(kpie,kpje,omask)
      &                      - refill*burial(i,j,issssil)
         burial(i,j,issster) = burial(i,j,issster)                      &
      &                      - refill*burial(i,j,issster)
-
+#ifdef cisonew
+        burial(i,j,issso13) = burial(i,j,issso13)                      &
+     &                      - refill*burial(i,j,issso13)
+        burial(i,j,isssc13) = burial(i,j,isssc13)                      &
+     &                      - refill*burial(i,j,isssc13)
+        burial(i,j,issso14) = burial(i,j,issso14)                      &
+     &                      - refill*burial(i,j,issso14)
+        burial(i,j,isssc14) = burial(i,j,isssc14)                      &
+     &                      - refill*burial(i,j,isssc14)
+#endif
       endif
       enddo !end i-loop
       enddo !end j-loop
@@ -262,7 +281,7 @@ SUBROUTINE SEDSHI(kpie,kpje,omask)
       enddo !end j-loop
 !$OMP END PARALLEL DO
 
-      do iv=1,4
+      do iv=1,nsedtra
 !$OMP PARALLEL DO PRIVATE(i,uebers,frac)
       do j=1,kpje
       do i=1,kpie
@@ -272,20 +291,6 @@ SUBROUTINE SEDSHI(kpie,kpje,omask)
           frac=porsol(i,j,k)*seddw(k)/(porsol(i,j,k-1)*seddw(k-1))
           sedlay(i,j,k,iv)=sedlay(i,j,k,iv)-uebers
           sedlay(i,j,k-1,iv)=sedlay(i,j,k-1,iv)+uebers*frac
-#ifdef cisonew
-          if(iv.eq.issso12)then
-            sedlay(i,j,k,issso13)  =sedlay(i,j,k,issso13)-uebers
-            sedlay(i,j,k-1,issso13)=sedlay(i,j,k-1,issso13)+uebers*frac
-            sedlay(i,j,k,issso14)  =sedlay(i,j,k,issso14)-uebers
-            sedlay(i,j,k-1,issso14)=sedlay(i,j,k-1,issso14)+uebers*frac
-          endif
-          if(iv.eq.isssc12)then
-            sedlay(i,j,k,isssc13)  =sedlay(i,j,k,isssc13)-uebers
-            sedlay(i,j,k-1,isssc13)=sedlay(i,j,k-1,isssc13)+uebers*frac
-            sedlay(i,j,k,isssc14)  =sedlay(i,j,k,isssc14)-uebers
-            sedlay(i,j,k-1,isssc14)=sedlay(i,j,k-1,isssc14)+uebers*frac
-          endif
-#endif
         endif
       enddo !end i-loop
       enddo !end j-loop