C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_ini_atmos.F,v 1.3 2015/12/25 00:48:45 jmc Exp $ C $Name: $ #include "DIC_OPTIONS.h" #include "PTRACERS_OPTIONS.h" CBOP C !ROUTINE: DIC_INI_ATMOS C !INTERFACE: ========================================================== SUBROUTINE DIC_INI_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_carbon(nSx,nSy) _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 C variables for reading CO2 input files _RL aWght, bWght _RL atm_pCO2 CEOP C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| ioUnit = standardMessageUnit C- initialise <-- done earlier, in dic_init_varia.F c total_atmos_carbon = 0. c atpco2 = 0. C user specified value (or default = 278 ppm)- set only once IF ( dic_int1.EQ.0 .OR. dic_int1.EQ.1 ) THEN atm_pCO2 = dic_pCO2 ELSEIF (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_carbon(bi,bj) = 0. 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_carbon, total_carbon, myThid ) C use dic_pCO2 as initial atmospheric pCO2 (not restart case): _BEGIN_MASTER(myThid) atpco2 = dic_pCO2 total_atmos_carbon = total_atmos_moles*dic_pCO2 _END_MASTER(myThid) C restart case: read previous atmospheric CO2 content & pCO2 from pickup file IF ( nIter0.GT.PTRACERS_Iter0 .OR. & (nIter0.EQ.PTRACERS_Iter0 .AND. pickupSuff.NE.' ') & ) THEN CALL DIC_READ_CO2_PICKUP( nIter0, myThid ) ENDIF _BEGIN_MASTER(myThid) C save initial content in common block var total_ocean_carbon = total_carbon atpco2 = total_atmos_carbon/total_atmos_moles C store initial content: total_ocean_carbon_start = total_carbon total_atmos_carbon_start = total_atmos_carbon total_ocean_carbon_old = total_carbon total_atmos_carbon_old = total_atmos_carbon total_ocean_carbon_year = total_carbon total_atmos_carbon_year = total_atmos_carbon C print out budget: WRITE(ioUnit,*) 'QQ atmos C, total, pCo2', & total_atmos_carbon, atpco2 total_carbon = total_atmos_carbon + total_ocean_carbon total_carbon_old = total_carbon carbon_diff = 0. WRITE(ioUnit,*) 'QQ total C, current, old, diff', & total_carbon, total_carbon_old, carbon_diff carbon_diff = 0. 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', & carbon_diff, carbon_diff _END_MASTER(myThid) _BARRIER atm_pCO2 = atpco2 ELSE atm_pCO2 = dic_pCO2 C end if dic_int1 = 0,1,2,or 3 ENDIF C-- Set initial AtmospCO2: 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 #endif /* ALLOW_DIC */ RETURN END