C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_atmos.F,v 1.20 2016/01/11 21:46:55 jmc Exp $
C $Name: $
#include "DIC_OPTIONS.h"
#include "PTRACERS_OPTIONS.h"
CBOP
C !ROUTINE: DIC_ATMOS
C !INTERFACE: ==========================================================
SUBROUTINE DIC_ATMOS( myTime, myIter, myThid )
C !DESCRIPTION:
C Calculate the atmospheric pCO2
C dic_int1:
C 0=use default 278.d-6
C 1=use constant value - dic_pCO2, read in from data.dic
C 2=read in from file
C 3=interact with atmospheric box (use dic_pCO2 as initial atmos. value)
C !USES: ===============================================================
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DIC_VARS.h"
#include "PTRACERS_SIZE.h"
#include "PTRACERS_PARAMS.h"
#include "PTRACERS_FIELDS.h"
#include "DIC_ATMOS.h"
C !INPUT PARAMETERS: ===================================================
C myTime :: current time
C myIter :: current iteration number
C myThid :: my Thread Id number
_RL myTime
INTEGER myIter, myThid
#ifdef ALLOW_DIC
C !FUNCTIONS: ====================================================
LOGICAL DIFFERENT_MULTIPLE
EXTERNAL
C !LOCAL VARIABLES: ====================================================
C total_atmos_moles :: atmosphere total gas content (should be parameter)
_RL total_atmos_moles
INTEGER bi, bj, i,j,k
INTEGER ntim
_RL tile_flux (nSx,nSy)
_RL tile_carbon(nSx,nSy)
_RL total_flux
_RL total_carbon
C for carbon budget ouput
INTEGER ioUnit
_RL total_ocean_carbon_old
_RL total_atmos_carbon_old
_RL total_carbon_old, carbon_diff
_RL year_diff_ocean, year_diff_atmos, year_total
_RL start_diff_ocean, start_diff_atmos, start_total
C variables for reading CO2 input files
_RL aWght, bWght
_RL atm_pCO2
LOGICAL timeCO2budget
CEOP
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ioUnit = standardMessageUnit
IF ( dic_int1.EQ.2 ) THEN
C read from a file and linearly interpolate between file entries
C (note: dic_int2=number entries to read
C dic_int3=start timestep,
C dic_int4=timestep between file entries)
ntim=int((myIter-dic_int3)/dic_int4)+1
aWght = FLOAT(myIter-dic_int3)
bWght = FLOAT(dic_int4)
aWght = 0.5 _d 0 + aWght/bWght - FLOAT(ntim-1)
IF (aWght.GT.1. _d 0) THEN
ntim=ntim+1
aWght=aWght-1. _d 0
ENDIF
bWght = 1. _d 0 - aWght
atm_pCO2 = co2atmos(ntim)*bWght + co2atmos(ntim+1)*aWght
WRITE(ioUnit,*) 'weights',ntim, aWght, bWght, atm_pCO2
ELSEIF ( dic_int1.EQ.3 ) THEN
C interactive atmosphere
C Mass dry atmosphere = (5.1352+/-0.0003)d18 kg (Trenberth & Smith,
C Journal of Climate 2005)
C and Mean molecular mass air = 28.97 g/mol (NASA earth fact sheet)
total_atmos_moles = 1.77 _d 20
C for 278ppmv we need total_atmos_carbon=4.9206e+16
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
tile_flux(bi,bj) = 0.
tile_carbon(bi,bj) = 0.
DO j=1,sNy
DO i=1,sNx
tile_flux(bi,bj) = tile_flux(bi,bj)
& + fluxCO2(i,j,bi,bj)*rA(i,j,bi,bj)
& *maskC(i,j,1,bi,bj)*dTtracerLev(1)
ENDDO
ENDDO
DO k=1,Nr
DO j=1,sNy
DO i=1,sNx
tile_carbon(bi,bj) = tile_carbon(bi,bj)
& + ( pTracer(i,j,k,bi,bj,1)
#ifdef DIC_BIOTIC
& +R_cp*pTracer(i,j,k,bi,bj,4)
#endif
& ) * rA(i,j,bi,bj)
& *drF(k)*hFacC(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
CALL GLOBAL_SUM_TILE_RL( tile_flux, total_flux, myThid )
CALL GLOBAL_SUM_TILE_RL( tile_carbon, total_carbon, myThid )
#ifdef ALLOW_AUTODIFF
C--NOTE: we DO need to make this section Master-Thread only to prevent multiple
C update (by each thread) of shared variable "total_atmos_carbon"
C- However, TAF does not recognize MITgcm multi-threading and shared variable
C status and see a conditional (if myThid=1) reset of "atpco2" that results
C in major recomputations.
#else
_BEGIN_MASTER(myThid)
#endif
C store previous content:
total_ocean_carbon_old = total_ocean_carbon
total_atmos_carbon_old = total_atmos_carbon
C calculate new atmos pCO2
total_atmos_carbon = total_atmos_carbon - total_flux
total_ocean_carbon = total_carbon
atpco2 = total_atmos_carbon/total_atmos_moles
WRITE(ioUnit,*) 'QQ atmos C, total, pCo2',
& total_atmos_carbon, atpco2
total_carbon = total_atmos_carbon + total_ocean_carbon
total_carbon_old = total_atmos_carbon_old+total_ocean_carbon_old
carbon_diff = total_carbon-total_carbon_old
WRITE(ioUnit,*) 'QQ total C, current, old, diff',
& total_carbon, total_carbon_old, carbon_diff
carbon_diff = total_ocean_carbon-total_ocean_carbon_old
WRITE(ioUnit,*) 'QQ ocean C, current, old, diff',
& total_ocean_carbon, total_ocean_carbon_old, carbon_diff
WRITE(ioUnit,*) 'QQ air-sea flux, addition diff',
& total_flux, carbon_diff-total_flux
C if end of forcing cycle, find total change in ocean carbon
timeCO2budget =
& DIFFERENT_MULTIPLE(externForcingCycle,myTime,deltaTClock)
IF ( timeCO2budget ) THEN
year_diff_ocean = total_ocean_carbon-total_ocean_carbon_year
year_diff_atmos = total_atmos_carbon-total_atmos_carbon_year
year_total = (total_ocean_carbon+total_atmos_carbon) -
& (total_ocean_carbon_year+total_atmos_carbon_year)
start_diff_ocean = total_ocean_carbon-total_ocean_carbon_start
start_diff_atmos = total_atmos_carbon-total_atmos_carbon_start
start_total = (total_ocean_carbon+total_atmos_carbon) -
& (total_ocean_carbon_start+total_atmos_carbon_start)
WRITE(ioUnit,*) 'QQ YEAR END'
WRITE(ioUnit,*) 'year diff: ocean, atmos, total',
& year_diff_ocean, year_diff_atmos, year_total
WRITE(ioUnit,*) 'start diff: ocean, atmos, total ',
& start_diff_ocean, start_diff_atmos, start_total
total_ocean_carbon_year = total_ocean_carbon
total_atmos_carbon_year = total_atmos_carbon
ENDIF
#ifndef ALLOW_AUTODIFF
_END_MASTER(myThid)
_BARRIER
#endif
atm_pCO2 = atpco2
ELSE
atm_pCO2 = dic_pCO2
ENDIF
#ifndef ALLOW_AUTODIFF
IF ( dic_int1.EQ.2 .OR. dic_int1.EQ.3 ) THEN
#endif
C-- Set AtmospCO2 for next iteration:
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
AtmospCO2(i,j,bi,bj) = atm_pCO2
ENDDO
ENDDO
ENDDO
ENDDO
#ifndef ALLOW_AUTODIFF
ENDIF
#endif
#endif /* ALLOW_DIC */
RETURN
END