pkg/dic/dic_ini_atmos.F

3f33260b4d Jean*0001 #include "DIC_OPTIONS.h"
                0002 #include "PTRACERS_OPTIONS.h"
                0003
                0004 CBOP
                0005 C !ROUTINE: DIC_INI_ATMOS
                0006
                0007 C !INTERFACE: ==========================================================
                0008       SUBROUTINE DIC_INI_ATMOS( myTime, myIter, myThid )
                0009
                0010 C !DESCRIPTION:
                0011 C  Calculate the atmospheric pCO2
                0012 C  dic_int1:
                0013 C  0=use default 278.d-6
                0014 C  1=use constant value - dic_pCO2, read in from data.dic
                0015 C  2=read in from file
                0016 C  3=interact with atmospheric box (use dic_pCO2 as initial atmos. value)
                0017
                0018 C !USES: ===============================================================
                0019       IMPLICIT NONE
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
                0023 #include "GRID.h"
                0024 #include "DIC_VARS.h"
                0025 #include "PTRACERS_SIZE.h"
                0026 #include "PTRACERS_PARAMS.h"
                0027 #include "PTRACERS_FIELDS.h"
                0028 #include "DIC_ATMOS.h"
                0029
                0030 C !INPUT PARAMETERS: ===================================================
                0031 C  myTime             :: current time
                0032 C  myIter             :: current iteration number
                0033 C  myThid             :: my Thread Id number
                0034       _RL myTime
                0035       INTEGER myIter, myThid
                0036
                0037 #ifdef ALLOW_DIC
                0038
                0039 C !FUNCTIONS:       ====================================================
                0040       LOGICAL  DIFFERENT_MULTIPLE
                0041       EXTERNAL DIFFERENT_MULTIPLE
                0042
                0043 C !LOCAL VARIABLES: ====================================================
                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_carbon(nSx,nSy)
                0050       _RL total_carbon
                0051
                0052 C for carbon budget ouput
                0053       INTEGER ioUnit
                0054       _RL total_ocean_carbon_old
                0055       _RL total_atmos_carbon_old
                0056       _RL total_carbon_old, carbon_diff
                0057 C variables for reading CO2 input files
                0058       _RL aWght, bWght
                0059       _RL atm_pCO2
                0060 CEOP
                0061
                0062 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0063
                0064       ioUnit = standardMessageUnit
                0065
                0066 C-    initialise <-- done earlier, in dic_init_varia.F
                0067 c     total_atmos_carbon = 0.
                0068 c     atpco2 = 0.
                0069
                0070 C user specified value (or default = 278 ppm)- set only once
                0071       IF ( dic_int1.EQ.0 .OR. dic_int1.EQ.1 ) THEN
                0072
                0073          atm_pCO2 = dic_pCO2
                0074
                0075       ELSEIF (dic_int1.EQ.2) THEN
                0076 C read from a file and linearly interpolate between file entries
                0077 C     (note:  dic_int2=number entries to read
                0078 C             dic_int3=start timestep,
                0079 C             dic_int4=timestep between file entries)
                0080
                0081         ntim=int((myIter-dic_int3)/dic_int4)+1
                0082         aWght = FLOAT(myIter-dic_int3)
                0083         bWght = FLOAT(dic_int4)
                0084         aWght = 0.5 _d 0 + aWght/bWght - FLOAT(ntim-1)
                0085         IF (aWght.GT.1. _d 0) THEN
                0086           ntim=ntim+1
                0087           aWght=aWght-1. _d 0
                0088         ENDIF
                0089         bWght = 1. _d 0 - aWght
                0090         atm_pCO2 = co2atmos(ntim)*bWght + co2atmos(ntim+1)*aWght
                0091         WRITE(ioUnit,*) 'weights',ntim, aWght, bWght, atm_pCO2
                0092
                0093       ELSEIF (dic_int1.EQ.3) THEN
                0094 C interactive atmosphere
                0095
                0096 C Mass dry atmosphere = (5.1352+/-0.0003)d18 kg (Trenberth & Smith,
                0097 C Journal of Climate 2005)
                0098 C and Mean molecular mass air = 28.97 g/mol (NASA earth fact sheet)
                0099        total_atmos_moles = 1.77 _d 20
                0100 C for 278ppmv we need total_atmos_carbon=4.9206e+16
                0101
                0102        DO bj=myByLo(myThid),myByHi(myThid)
                0103         DO bi=myBxLo(myThid),myBxHi(myThid)
                0104          tile_carbon(bi,bj) = 0.
                0105          DO k=1,Nr
                0106           DO j=1,sNy
                0107            DO i=1,sNx
                0108              tile_carbon(bi,bj) = tile_carbon(bi,bj)
                0109      &            + ( pTracer(i,j,k,bi,bj,1)
                0110 #ifdef DIC_BIOTIC
                0111      &               +R_cp*pTracer(i,j,k,bi,bj,4)
                0112 #endif
                0113      &              ) * rA(i,j,bi,bj)
                0114      &                *drF(k)*hFacC(i,j,k,bi,bj)
                0115            ENDDO
                0116           ENDDO
                0117          ENDDO
                0118         ENDDO
                0119        ENDDO
                0120
                0121        CALL GLOBAL_SUM_TILE_RL( tile_carbon, total_carbon, myThid )
                0122
                0123 C use dic_pCO2 as initial atmospheric pCO2 (not restart case):
                0124          _BEGIN_MASTER(myThid)
                0125          atpco2 = dic_pCO2
                0126          total_atmos_carbon = total_atmos_moles*dic_pCO2
                0127          _END_MASTER(myThid)
                0128 C restart case: read previous atmospheric CO2 content & pCO2 from pickup file
                0129          IF ( nIter0.GT.PTRACERS_Iter0 .OR.
                0130      &       (nIter0.EQ.PTRACERS_Iter0 .AND. pickupSuff.NE.' ')
                0131      &      ) THEN
                0132            CALL DIC_READ_CO2_PICKUP( nIter0, myThid )
                0133          ENDIF
                0134
                0135        _BEGIN_MASTER(myThid)
                0136 C save initial content in common block var
                0137         total_ocean_carbon = total_carbon
                0138         atpco2 = total_atmos_carbon/total_atmos_moles
                0139
                0140 C store initial content:
                0141         total_ocean_carbon_start = total_carbon
                0142         total_atmos_carbon_start = total_atmos_carbon
                0143         total_ocean_carbon_old   = total_carbon
                0144         total_atmos_carbon_old   = total_atmos_carbon
                0145         total_ocean_carbon_year  = total_carbon
                0146         total_atmos_carbon_year  = total_atmos_carbon
                0147
                0148 C print out budget:
                0149         WRITE(ioUnit,*) 'QQ atmos C, total, pCo2',
                0150      &                     total_atmos_carbon, atpco2
                0151         total_carbon = total_atmos_carbon + total_ocean_carbon
                0152         total_carbon_old = total_carbon
                0153         carbon_diff = 0.
                0154         WRITE(ioUnit,*) 'QQ total C, current, old, diff',
                0155      &                     total_carbon, total_carbon_old, carbon_diff
                0156         carbon_diff = 0.
                0157         WRITE(ioUnit,*) 'QQ ocean C, current, old, diff',
                0158      &         total_ocean_carbon, total_ocean_carbon_old, carbon_diff
                0159         WRITE(ioUnit,*) 'QQ air-sea flux, addition diff',
64d6173181 Jean*0160      &                     carbon_diff, carbon_diff
3f33260b4d Jean*0161
                0162        _END_MASTER(myThid)
                0163        _BARRIER
                0164
                0165         atm_pCO2 = atpco2
                0166       ELSE
                0167         atm_pCO2 = dic_pCO2
                0168 C end if dic_int1 = 0,1,2,or 3
                0169       ENDIF
                0170
64d6173181 Jean*0171 C--   Set initial AtmospCO2:
3f33260b4d Jean*0172       DO bj=myByLo(myThid),myByHi(myThid)
                0173        DO bi=myBxLo(myThid),myBxHi(myThid)
                0174          DO j=1-OLy,sNy+OLy
                0175           DO i=1-OLx,sNx+OLx
                0176             AtmospCO2(i,j,bi,bj) = atm_pCO2
                0177           ENDDO
                0178          ENDDO
                0179        ENDDO
                0180       ENDDO
                0181
                0182 #endif /* ALLOW_DIC */
                0183
                0184       RETURN
                0185       END