pkg/dic/dic_atmos.F

1ef1c61fb1 Step*0001 #include "DIC_OPTIONS.h"
                0002 #include "PTRACERS_OPTIONS.h"
                0003
                0004 CBOP
                0005 C !ROUTINE: DIC_ATMOS
                0006
                0007 C !INTERFACE: ==========================================================
3f33260b4d Jean*0008       SUBROUTINE DIC_ATMOS( myTime, myIter, myThid )
1ef1c61fb1 Step*0009
                0010 C !DESCRIPTION:
                0011 C  Calculate the atmospheric pCO2
254eed5bf2 Davi*0012 C  dic_int1:
1ef1c61fb1 Step*0013 C  0=use default 278.d-6
254eed5bf2 Davi*0014 C  1=use constant value - dic_pCO2, read in from data.dic
f4d7963907 Jean*0015 C  2=read in from file
175a18b00a Jean*0016 C  3=interact with atmospheric box (use dic_pCO2 as initial atmos. value)
                0017
1ef1c61fb1 Step*0018 C !USES: ===============================================================
                0019       IMPLICIT NONE
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
                0023 #include "GRID.h"
2ef8966791 Davi*0024 #include "DIC_VARS.h"
df43453050 Davi*0025 #include "PTRACERS_SIZE.h"
175a18b00a Jean*0026 #include "PTRACERS_PARAMS.h"
df43453050 Davi*0027 #include "PTRACERS_FIELDS.h"
1ef1c61fb1 Step*0028 #include "DIC_ATMOS.h"
                0029
                0030 C !INPUT PARAMETERS: ===================================================
3f33260b4d Jean*0031 C  myTime             :: current time
                0032 C  myIter             :: current iteration number
                0033 C  myThid             :: my Thread Id number
1ef1c61fb1 Step*0034       _RL myTime
175a18b00a Jean*0035       INTEGER myIter, myThid
1ef1c61fb1 Step*0036
175a18b00a Jean*0037 #ifdef ALLOW_DIC
                0038
                0039 C !FUNCTIONS:       ====================================================
1ef1c61fb1 Step*0040       LOGICAL  DIFFERENT_MULTIPLE
                0041       EXTERNAL DIFFERENT_MULTIPLE
                0042
                0043 C !LOCAL VARIABLES: ====================================================
175a18b00a Jean*0044 C   total_atmos_moles :: atmosphere total gas content (should be parameter)
                0045       _RL total_atmos_moles
                0046       INTEGER bi, bj, i,j,k
                0047       INTEGER ntim
                0048
                0049       _RL tile_flux  (nSx,nSy)
                0050       _RL tile_carbon(nSx,nSy)
                0051       _RL total_flux
                0052       _RL total_carbon
                0053
                0054 C for carbon budget ouput
                0055       INTEGER ioUnit
                0056       _RL total_ocean_carbon_old
                0057       _RL total_atmos_carbon_old
                0058       _RL total_carbon_old, carbon_diff
                0059       _RL year_diff_ocean, year_diff_atmos, year_total
                0060       _RL start_diff_ocean, start_diff_atmos, start_total
1ef1c61fb1 Step*0061 C variables for reading CO2 input files
                0062       _RL aWght, bWght
3f33260b4d Jean*0063       _RL atm_pCO2
175a18b00a Jean*0064       LOGICAL timeCO2budget
                0065 CEOP
1ef1c61fb1 Step*0066
175a18b00a Jean*0067 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1ef1c61fb1 Step*0068
175a18b00a Jean*0069       ioUnit = standardMessageUnit
1ef1c61fb1 Step*0070
3f33260b4d Jean*0071       IF ( dic_int1.EQ.2 ) THEN
                0072 C read from a file and linearly interpolate between file entries
                0073 C     (note:  dic_int2=number entries to read
                0074 C             dic_int3=start timestep,
                0075 C             dic_int4=timestep between file entries)
175a18b00a Jean*0076         ntim=int((myIter-dic_int3)/dic_int4)+1
                0077         aWght = FLOAT(myIter-dic_int3)
                0078         bWght = FLOAT(dic_int4)
                0079         aWght = 0.5 _d 0 + aWght/bWght - FLOAT(ntim-1)
                0080         IF (aWght.GT.1. _d 0) THEN
                0081           ntim=ntim+1
                0082           aWght=aWght-1. _d 0
                0083         ENDIF
                0084         bWght = 1. _d 0 - aWght
3f33260b4d Jean*0085         atm_pCO2 = co2atmos(ntim)*bWght + co2atmos(ntim+1)*aWght
                0086         WRITE(ioUnit,*) 'weights',ntim, aWght, bWght, atm_pCO2
1ef1c61fb1 Step*0087
3f33260b4d Jean*0088       ELSEIF ( dic_int1.EQ.3 ) THEN
175a18b00a Jean*0089 C interactive atmosphere
1ef1c61fb1 Step*0090
175a18b00a Jean*0091 C Mass dry atmosphere = (5.1352+/-0.0003)d18 kg (Trenberth & Smith,
                0092 C Journal of Climate 2005)
                0093 C and Mean molecular mass air = 28.97 g/mol (NASA earth fact sheet)
3f33260b4d Jean*0094        total_atmos_moles = 1.77 _d 20
175a18b00a Jean*0095 C for 278ppmv we need total_atmos_carbon=4.9206e+16
1ef1c61fb1 Step*0096
                0097        DO bj=myByLo(myThid),myByHi(myThid)
                0098         DO bi=myBxLo(myThid),myBxHi(myThid)
175a18b00a Jean*0099          tile_flux(bi,bj)   = 0.
                0100          tile_carbon(bi,bj) = 0.
60a43582f7 Jean*0101          DO j=1,sNy
                0102            DO i=1,sNx
                0103              tile_flux(bi,bj) = tile_flux(bi,bj)
                0104      &                        + fluxCO2(i,j,bi,bj)*rA(i,j,bi,bj)
                0105      &                         *maskC(i,j,1,bi,bj)*dTtracerLev(1)
                0106            ENDDO
                0107          ENDDO
175a18b00a Jean*0108          DO k=1,Nr
1ef1c61fb1 Step*0109           DO j=1,sNy
175a18b00a Jean*0110            DO i=1,sNx
                0111              tile_carbon(bi,bj) = tile_carbon(bi,bj)
                0112      &            + ( pTracer(i,j,k,bi,bj,1)
df43453050 Davi*0113 #ifdef DIC_BIOTIC
175a18b00a Jean*0114      &               +R_cp*pTracer(i,j,k,bi,bj,4)
df43453050 Davi*0115 #endif
175a18b00a Jean*0116      &              ) * rA(i,j,bi,bj)
                0117      &                *drF(k)*hFacC(i,j,k,bi,bj)
                0118            ENDDO
1ef1c61fb1 Step*0119           ENDDO
                0120          ENDDO
                0121         ENDDO
                0122        ENDDO
e5151264e4 Step*0123
175a18b00a Jean*0124        CALL GLOBAL_SUM_TILE_RL( tile_flux,   total_flux,   myThid )
                0125        CALL GLOBAL_SUM_TILE_RL( tile_carbon, total_carbon, myThid )
                0126
3f33260b4d Jean*0127 #ifdef ALLOW_AUTODIFF
                0128 C--NOTE: we DO need to make this section Master-Thread only to prevent multiple
                0129 C   update (by each thread) of shared variable "total_atmos_carbon"
                0130 C- However, TAF does not recognize MITgcm multi-threading and shared variable
                0131 C   status and see a conditional (if myThid=1) reset of "atpco2" that results
                0132 C   in major recomputations.
                0133 #else
                0134        _BEGIN_MASTER(myThid)
79e00a2386 Step*0135 #endif
175a18b00a Jean*0136 C store previous content:
3f33260b4d Jean*0137        total_ocean_carbon_old = total_ocean_carbon
                0138        total_atmos_carbon_old = total_atmos_carbon
175a18b00a Jean*0139 C calculate new atmos pCO2
3f33260b4d Jean*0140        total_atmos_carbon = total_atmos_carbon - total_flux
175a18b00a Jean*0141        total_ocean_carbon = total_carbon
                0142        atpco2 = total_atmos_carbon/total_atmos_moles
1ef1c61fb1 Step*0143
175a18b00a Jean*0144        WRITE(ioUnit,*) 'QQ atmos C, total, pCo2',
                0145      &                     total_atmos_carbon, atpco2
3f33260b4d Jean*0146        total_carbon = total_atmos_carbon + total_ocean_carbon
                0147        total_carbon_old = total_atmos_carbon_old+total_ocean_carbon_old
                0148        carbon_diff = total_carbon-total_carbon_old
175a18b00a Jean*0149        WRITE(ioUnit,*) 'QQ total C, current, old, diff',
                0150      &                     total_carbon, total_carbon_old, carbon_diff
3f33260b4d Jean*0151        carbon_diff = total_ocean_carbon-total_ocean_carbon_old
175a18b00a Jean*0152        WRITE(ioUnit,*) 'QQ ocean C, current, old, diff',
                0153      &         total_ocean_carbon, total_ocean_carbon_old, carbon_diff
                0154        WRITE(ioUnit,*) 'QQ air-sea flux, addition diff',
                0155      &                     total_flux, carbon_diff-total_flux
                0156
                0157 C if end of forcing cycle, find total change in ocean carbon
3f33260b4d Jean*0158        timeCO2budget =
                0159      &    DIFFERENT_MULTIPLE(externForcingCycle,myTime,deltaTClock)
                0160        IF ( timeCO2budget ) THEN
175a18b00a Jean*0161           year_diff_ocean = total_ocean_carbon-total_ocean_carbon_year
                0162           year_diff_atmos = total_atmos_carbon-total_atmos_carbon_year
                0163           year_total = (total_ocean_carbon+total_atmos_carbon) -
                0164      &               (total_ocean_carbon_year+total_atmos_carbon_year)
                0165           start_diff_ocean = total_ocean_carbon-total_ocean_carbon_start
                0166           start_diff_atmos = total_atmos_carbon-total_atmos_carbon_start
                0167           start_total = (total_ocean_carbon+total_atmos_carbon) -
                0168      &               (total_ocean_carbon_start+total_atmos_carbon_start)
                0169           WRITE(ioUnit,*) 'QQ YEAR END'
                0170           WRITE(ioUnit,*) 'year diff: ocean, atmos, total',
                0171      &                  year_diff_ocean,  year_diff_atmos,  year_total
                0172           WRITE(ioUnit,*) 'start diff: ocean, atmos, total ',
                0173      &                 start_diff_ocean, start_diff_atmos, start_total
                0174
                0175           total_ocean_carbon_year = total_ocean_carbon
                0176           total_atmos_carbon_year = total_atmos_carbon
                0177        ENDIF
                0178
3f33260b4d Jean*0179 #ifndef ALLOW_AUTODIFF
175a18b00a Jean*0180        _END_MASTER(myThid)
                0181        _BARRIER
3f33260b4d Jean*0182 #endif
                0183
                0184        atm_pCO2 = atpco2
                0185       ELSE
                0186        atm_pCO2 = dic_pCO2
                0187       ENDIF
175a18b00a Jean*0188
3f33260b4d Jean*0189 #ifndef ALLOW_AUTODIFF
6569fc5765 Jean*0190       IF ( dic_int1.EQ.2 .OR. dic_int1.EQ.3 ) THEN
3f33260b4d Jean*0191 #endif
175a18b00a Jean*0192 C--    Set AtmospCO2 for next iteration:
1ef1c61fb1 Step*0193        DO bj=myByLo(myThid),myByHi(myThid)
                0194         DO bi=myBxLo(myThid),myBxHi(myThid)
                0195          DO j=1-OLy,sNy+OLy
                0196           DO i=1-OLx,sNx+OLx
3f33260b4d Jean*0197             AtmospCO2(i,j,bi,bj) = atm_pCO2
1ef1c61fb1 Step*0198           ENDDO
                0199          ENDDO
                0200         ENDDO
caaf66502e Davi*0201        ENDDO
3f33260b4d Jean*0202 #ifndef ALLOW_AUTODIFF
175a18b00a Jean*0203       ENDIF
3f33260b4d Jean*0204 #endif
1ef1c61fb1 Step*0205
175a18b00a Jean*0206 #endif /* ALLOW_DIC */
1ef1c61fb1 Step*0207
5f38feeebb Jean*0208       RETURN
                0209       END