C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_biotic_forcing.F,v 1.32 2016/01/16 21:57:12 jmc Exp $ C $Name: $ #include "DIC_OPTIONS.h" CBOP C !ROUTINE: DIC_BIOTIC_FORCING C !INTERFACE: ========================================================== SUBROUTINE DIC_BIOTIC_FORCING( U PTR_DIC, PTR_ALK, PTR_PO4, PTR_DOP, #ifdef ALLOW_O2 U PTR_O2, #endif #ifdef ALLOW_FE U PTR_FE, #endif I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) C !DESCRIPTION: C updates all the tracers for the effects of air-sea exchange, biological c activity and remineralization C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "DYNVARS.h" #include "DIC_VARS.h" #include "PTRACERS_SIZE.h" #include "PTRACERS_PARAMS.h" C !INPUT/OUTPUT PARAMETERS: =================================================== C PTR_DIC :: dissolced inorganic carbon C PTR_ALK :: alkalinity C PTR_PO4 :: phosphate c PTR_DOP :: dissolve organic phosphurous c PTR_O2 :: oxygen C PTR_FE :: iron c bi, bj :: current tile indices C myIter :: current timestep C myTime :: current time C myThid :: thread number _RL PTR_DIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL PTR_ALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL PTR_DOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) #ifdef ALLOW_O2 _RL PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) #endif #ifdef ALLOW_FE _RL PTR_FE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) #endif INTEGER bi, bj, iMin, iMax, jMin, jMax INTEGER myIter _RL myTime INTEGER myThid #ifdef ALLOW_PTRACERS #ifdef DIC_BIOTIC C !LOCAL VARIABLES: ==================================================== C i,j,k :: loop indices C G* :: tendency term for the tracers C SURA :: tendency of alkalinity due to freshwater C SURC :: tendency of DIC due to air-sea exchange C and virtual flux C SURO :: tendency of O2 due to air-sea exchange C GPO4 :: tendency of PO4 due to biological productivity, C exchange with DOP pool and reminerization C CAR :: carbonate changes due to biological C productivity and remineralization C BIOac :: biological productivity C RDOP :: DOP sink due to remineralization C pflux :: changes to PO4 due to flux and remineralization C CAR_S :: carbonate sink C cflux :: carbonate changes due to flux and remineralization C freefe :: iron not bound to ligand _RL GDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL GALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL GPO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL GDOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL SURA(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL SURC(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL SURO(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL CAR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL BIOac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL RDOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL pflux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL exportflux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL CAR_S(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL cflux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) #ifdef ALLOW_O2 _RL GO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) #endif #ifdef ALLOW_FE _RL GFE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL freefe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) #endif INTEGER i,j,k #ifdef CAR_DISS INTEGER nCALCITEstep #endif #ifdef ALLOW_FE # ifdef SEDFE INTEGER kBottom # endif #endif CEOP #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_ENTER('DIC_BIOTIC_FORCING',myThid) #endif IF ( useThSIce .OR. useSEAICE .OR. useCoupler ) THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('DIC_FIELDS_UPDATE',myThid) #endif CALL DIC_FIELDS_UPDATE( I bi, bj, myTime, myIter, myThid ) ENDIF DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx RDOP(i,j,k) =0. _d 0 GDIC(i,j,k) =0. _d 0 GALK(i,j,k) =0. _d 0 GPO4(i,j,k) =0. _d 0 GDOP(i,j,k) =0. _d 0 CAR(i,j,k) =0. _d 0 BIOac(i,j,k) =0. _d 0 pflux(i,j,k) =0. _d 0 exportflux(i,j,k)=0. _d 0 cflux(i,j,k) =0. _d 0 CAR_S(i,j,k) =0. _d 0 #ifdef ALLOW_O2 GO2(i,j,k) =0. _d 0 #endif #ifdef ALLOW_FE GFE(i,j,k) =0. _d 0 C no longer needed after adding full initialisation of freefe in S/R FE_CHEM c freefe(i,j,k) =0. _d 0 #endif ENDDO ENDDO ENDDO DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx SURA(i,j) =0. _d 0 SURC(i,j) =0. _d 0 SURO(i,j) =0. _d 0 ENDDO ENDDO C carbon air-sea interaction #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('DIC_SURFFORCING',myThid) #endif CALL DIC_SURFFORCING( I PTR_DIC, PTR_ALK, PTR_PO4, O SURC, I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) C alkalinity air-sea interaction #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('ALK_SURFFORCING',myThid) #endif CALL ALK_SURFFORCING( I PTR_ALK, O SURA, I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) #ifdef ALLOW_O2 C oxygen air-sea interaction #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('O2_SURFFORCING',myThid) #endif CALL O2_SURFFORCING( I PTR_O2, O SURO, I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) #endif #ifdef ALLOW_FE C find free iron #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('FE_CHEM',myThid) #endif CALL FE_CHEM( bi, bj, iMin, iMax, jMin, jMax, U PTR_FE, O freefe, I myIter, myThid ) #endif C biological activity #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('BIO_EXPORT',myThid) #endif CALL BIO_EXPORT( I PTR_PO4, #ifdef ALLOW_FE I PTR_FE, #endif O BIOac, I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) C flux of po4 from layers with biological activity #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('PHOS_FLUX',myThid) #endif CALL PHOS_FLUX( I BIOac, U pflux, exportflux, I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) C- Carbonate sink DO k=1,Nr DO j=jMin,jMax DO i=iMin,iMax CAR_S(i,j,k)=BIOac(i,j,k)*R_CP*rain_ratio(i,j,bi,bj)* & (1. _d 0-DOPfraction) ENDDO ENDDO ENDDO C carbonate #ifdef CAR_DISS C dissolution only below saturation horizon C code following method by Karsten Friis nCALCITEstep = 3600 IF(myIter .lt. (nIter0+5) .or. & mod(myIter,nCALCITEstep) .eq. 0)THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('CALCITE_SATURATION',myThid) #endif CALL CALCITE_SATURATION( I PTR_DIC, PTR_ALK, PTR_PO4, I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) ENDIF #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('CAR_FLUX_OMEGA_TOP',myThid) #endif CALL CAR_FLUX_OMEGA_TOP( I BIOac, O cflux, I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) #else C old OCMIP way #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('CAR_FLUX',myThid) #endif CALL CAR_FLUX( I CAR_S, U cflux, I bi, bj, iMin, iMax, jMin, jMax, I myIter, myTime, myThid ) #endif C add all tendencies for PO4, DOP, ALK, DIC DO k=1,Nr DO j=jMin,jMax DO i=iMin,iMax #ifdef DIC_NO_NEG RDOP(i,j,k)= MAX(maskC(i,j,k,bi,bj)*KDOPRemin*PTR_DOP(i,j,k) & ,0. _d 0) #else RDOP(i,j,k)= maskC(i,j,k,bi,bj)*KDOPRemin*PTR_DOP(i,j,k) #endif GPO4(i,j,k)=-BIOac(i,j,k)+pflux(i,j,k) + RDOP(i,j,k) car(i,j,k) = cflux(i,j,k) - CAR_S(i,j,k) GDOP(i,j,k)=+BIOac(i,j,k)*DOPfraction - RDOP(i,j,k) GALK(i,j,k)=+2. _d 0 *car(i,j,k)-R_NP*GPO4(i,j,k) GDIC(i,j,k)=car(i,j,k)+R_CP*GPO4(i,j,k) #ifdef ALLOW_O2 if (PTR_O2(i,j,k).GT.O2crit) then GO2(i,j,k)= R_OP*GPO4(i,j,k) else GO2(i,j,k)= 0. _d 0 endif #endif #ifdef ALLOW_FE GFE(i,j,k) = R_FeP*GPO4(i,j,k) & -Kscav*freefe(i,j,k) #endif ENDDO ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax GALK(i,j,1)=GALK(i,j,1)+SURA(i,j) GDIC(i,j,1)=GDIC(i,j,1)+SURC(i,j) #ifdef ALLOW_O2 GO2(i,j,1) =GO2(i,j,1)+SURO(i,j) #endif #ifdef ALLOW_FE GFE(i,j,1)=GFE(i,j,1)+alpfe* & InputFe(i,j,bi,bj)*recip_drF(1) & *recip_hFacC(i,j,1,bi,bj) # ifdef SEDFE C include iron sediment source using the flux of po4 into bottom layer kBottom = MAX(kLowC(i,j,bi,bj),1) GFE(i,j,kBottom)=GFE(i,j,kBottom) & +( fesedflux_pcm*pflux(i,j,kBottom) + FeIntSec ) & *recip_drF(kBottom)*recip_hFacC(i,j,kBottom,bi,bj) # endif #endif ENDDO ENDDO IF ( useOBCS ) THEN DO k=1,Nr DO j=jMin,jMax DO i=iMin,iMax GDIC(i,j,k) = GDIC(i,j,k)*maskInC(i,j,bi,bj) GALK(i,j,k) = GALK(i,j,k)*maskInC(i,j,bi,bj) GPO4(i,j,k) = GPO4(i,j,k)*maskInC(i,j,bi,bj) GDOP(i,j,k) = GDOP(i,j,k)*maskInC(i,j,bi,bj) #ifdef ALLOW_O2 GO2(i,j,k) = GO2(i,j,k)*maskInC(i,j,bi,bj) #endif #ifdef ALLOW_FE GFE(i,j,k) = GFE(i,j,k)*maskInC(i,j,bi,bj) #endif ENDDO ENDDO ENDDO ENDIF C update DO k=1,Nr DO j=jMin,jMax DO i=iMin,iMax PTR_DIC(i,j,k)= & PTR_DIC(i,j,k)+GDIC(i,j,k)*PTRACERS_dTLev(k) PTR_ALK(i,j,k)= & PTR_ALK(i,j,k)+GALK(i,j,k)*PTRACERS_dTLev(k) PTR_PO4(i,j,k)= & PTR_PO4(i,j,k)+GPO4(i,j,k)*PTRACERS_dTLev(k) PTR_DOP(i,j,k)= & PTR_DOP(i,j,k)+GDOP(i,j,k)*PTRACERS_dTLev(k) #ifdef ALLOW_O2 PTR_O2(i,j,k)= & PTR_O2(i,j,k)+GO2(i,j,k)*PTRACERS_dTLev(k) #endif #ifdef ALLOW_FE PTR_FE(i,j,k)= & PTR_FE(i,j,k)+GFE(i,j,k)*PTRACERS_dTLev(k) #endif ENDDO ENDDO ENDDO #ifdef ALLOW_FE #ifdef MINFE c find free iron and get rid of insoluble part #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('FE_CHEM',myThid) #endif CALL FE_CHEM( bi, bj, iMin, iMax, jMin, jMax, U PTR_FE, O freefe, I myIter, myThid ) #endif #endif #ifdef ALLOW_TIMEAVE C save averages IF ( PTRACERS_taveFreq.GT.0. ) THEN DO k=1,Nr DO j=jMin,jMax DO i=iMin,iMax BIOave(i,j,k,bi,bj) =BIOave(i,j,k,bi,bj)+ & BIOac(i,j,k)*deltaTClock CARave(i,j,k,bi,bj) =CARave(i,j,k,bi,bj)+ & CAR(i,j,k)*deltaTClock OmegaCave(i,j,k,bi,bj)=OmegaCave(i,j,k,bi,bj)+ & OmegaC(i,j,k,bi,bj)*deltaTClock pfluxave(i,j,k,bi,bj) =pfluxave(i,j,k,bi,bj) + & pflux(i,j,k)*deltaTClock epfluxave(i,j,k,bi,bj)=epfluxave(i,j,k,bi,bj) + & exportflux(i,j,k)*deltaTClock cfluxave(i,j,k,bi,bj) =cfluxave(i,j,k,bi,bj) + & cflux(i,j,k)*deltaTClock ENDDO ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax SURave(i,j,bi,bj) =SURave(i,j,bi,bj)+ & SURC(i,j)*deltaTClock #ifdef ALLOW_O2 SUROave(i,j,bi,bj) =SUROave(i,j,bi,bj)+ & SURO(i,j)*deltaTClock #endif pCO2ave(i,j,bi,bj) =pCO2ave(i,j,bi,bj)+ & pCO2(i,j,bi,bj)*deltaTClock pHave(i,j,bi,bj) =pHave(i,j,bi,bj)+ & pH(i,j,bi,bj)*deltaTClock fluxCO2ave(i,j,bi,bj)=fluxCO2ave(i,j,bi,bj)+ & fluxCO2(i,j,bi,bj)*deltaTClock ENDDO ENDDO DIC_timeAve(bi,bj) = DIC_timeAve(bi,bj)+deltaTClock ENDIF #endif /* ALLOW_TIMEAVE*/ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN CALL DIAGNOSTICS_FILL(BIOac ,'DICBIOA ',0,Nr,2,bi,bj,myThid) CALL DIAGNOSTICS_FILL(CAR ,'DICCARB ',0,Nr,2,bi,bj,myThid) CALL DIAGNOSTICS_FILL(pCO2 ,'DICPCO2 ',0,1 ,1,bi,bj,myThid) CALL DIAGNOSTICS_FILL(fluxCO2,'DICCFLX ',0,1 ,1,bi,bj,myThid) CALL DIAGNOSTICS_FILL(pH ,'DICPHAV ',0,1 ,1,bi,bj,myThid) CALL DIAGNOSTICS_FILL(SURC ,'DICTFLX ',0,1 ,2,bi,bj,myThid) #ifdef ALLOW_O2 CALL DIAGNOSTICS_FILL(SURO ,'DICOFLX ',0,1 ,2,bi,bj,myThid) #endif ENDIF #endif /* ALLOW_DIAGNOSTICS */ #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_LEAVE('DIC_BIOTIC_FORCING',myThid) #endif #endif /* DIC_BIOTIC */ #endif /* ALLOW_PTRACERS */ RETURN END