/[MITgcm]/MITgcm/pkg/dic/dic_atmos.F
ViewVC logotype

Diff of /MITgcm/pkg/dic/dic_atmos.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.9 by dfer, Mon Apr 7 20:31:16 2008 UTC revision 1.16 by jmc, Sun Apr 11 21:04:39 2010 UTC
# Line 12  C !INTERFACE: ========================== Line 12  C !INTERFACE: ==========================
12    
13  C !DESCRIPTION:  C !DESCRIPTION:
14  C  Calculate the atmospheric pCO2  C  Calculate the atmospheric pCO2
15  C  gchem_int1:  C  dic_int1:
16  C  0=use default 278.d-6  C  0=use default 278.d-6
17  C  1=use constant value - gchem_rl1, read in from data.gchem  C  1=use constant value - dic_pCO2, read in from data.dic
18  C  2=read in from file  C  2=read in from file
19  C  3=interact with atmospheric box  C  3=interact with atmospheric box (use dic_pCO2 as initial atmos. value)
20    
21  C !USES: ===============================================================  C !USES: ===============================================================
22        IMPLICIT NONE        IMPLICIT NONE
23  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26  #include "GRID.h"  #include "GRID.h"
 #include "FFIELDS.h"  
27  #include "DIC_VARS.h"  #include "DIC_VARS.h"
28  #include "PTRACERS_SIZE.h"  #include "PTRACERS_SIZE.h"
29    #include "PTRACERS_PARAMS.h"
30  #include "PTRACERS_FIELDS.h"  #include "PTRACERS_FIELDS.h"
 #include "GCHEM.h"  
31  #include "DIC_ATMOS.h"  #include "DIC_ATMOS.h"
32    
33  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
 C  myThid               :: thread number  
 C  myIter               :: current timestep  
 C  myTime               :: current time  
34  C  istate               :: 0=initial call, 1=subsequent calls  C  istate               :: 0=initial call, 1=subsequent calls
35        INTEGER myIter, myThid, istate  C  myTime               :: current time
36    C  myIter               :: current iteration number
37    C  myThid               :: my Thread Id number
38          INTEGER istate
39        _RL myTime        _RL myTime
40          INTEGER myIter, myThid
41    
42    #ifdef ALLOW_DIC
43    
44  #ifdef ALLOW_PTRACERS  #ifdef USE_ATMOSCO2
45    C if coupled to atmsopheric model, use the
46    C CO2 value passed from the coupler
47    
48    #else /* USE_ATMOSCO2 */
49    
50    C !FUNCTIONS:       ====================================================
51        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
52        EXTERNAL DIFFERENT_MULTIPLE        EXTERNAL DIFFERENT_MULTIPLE
53    
54  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
55         INTEGER bi, bj, I,J,k  C   total_atmos_moles :: atmosphere total gas content (should be parameter)
56         INTEGER it, ntim        _RL total_atmos_moles
57  c        INTEGER bi, bj, i,j,k
58         _RL total_flux        INTEGER ntim
59         _RL total_ocean_carbon_old  
60         _RL total_atmos_carbon_old        _RL tile_flux  (nSx,nSy)
61         _RL total_atmos_moles        _RL tile_carbon(nSx,nSy)
62         _RL atpco2        _RL total_flux
63         _RL total_carbon_old, total_carbon, carbon_diff        _RL total_carbon
64         _RL tmp  
65         _RL year_diff_ocean, year_diff_atmos, year_total  C for carbon budget ouput
66         _RL start_diff_ocean, start_diff_atmos, start_total        INTEGER ioUnit
67          _RL total_ocean_carbon_old
68          _RL total_atmos_carbon_old
69          _RL total_carbon_old, carbon_diff
70          _RL year_diff_ocean, year_diff_atmos, year_total
71          _RL start_diff_ocean, start_diff_atmos, start_total
72  C variables for reading CO2 input files  C variables for reading CO2 input files
73          _RL tmp
74        _RL aWght, bWght        _RL aWght, bWght
 c  
        CHARACTER*(MAX_LEN_FNAM) fn  
        LOGICAL permCheckPoint  
 CEOP  
   
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
75    
76  c if coupled to atmsopheric model, use the        LOGICAL timeCO2budget
77  c Co2 value passed from the coupler  CEOP
 #ifndef USE_ATMOSCO2  
   
       IF ( nThreads .GT. 1 .AND.  
      &     ( gchem_int1.EQ.2 .OR. gchem_int1.EQ.3 ) ) THEN  
 C     Problem with I/O and global-sum for multi-threaded execution  
 C     Needs to be fixed before using this S/R in multi-threaded run  
         STOP 'S/R DIC_ATMOS: multi-threaded not right'  
       ENDIF  
   
 c default - set only once  
       if (gchem_int1.eq.0.and.istate.eq.0) then  
        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)=278.0 _d -6  
           ENDDO  
          ENDDO  
78    
79           ENDDO  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
         ENDDO  
       endif  
80    
81  c user specified value - set only once        ioUnit = standardMessageUnit
       if (gchem_int1.eq.1.and.istate.eq.0) then  
        DO bj=myByLo(myThid),myByHi(myThid)  
         DO bi=myBxLo(myThid),myBxHi(myThid)  
82    
83           DO j=1-OLy,sNy+OLy  C user specified value (or default = 278 ppm)- set only once
84            DO i=1-OLx,sNx+OLx        IF ( (dic_int1.EQ.0 .OR. dic_int1.EQ.1) .AND. istate.EQ.0 ) THEN
85               AtmospCO2(i,j,bi,bj)=gchem_rl1          DO bj=myByLo(myThid),myByHi(myThid)
86             DO bi=myBxLo(myThid),myBxHi(myThid)
87              DO j=1-OLy,sNy+OLy
88               DO i=1-OLx,sNx+OLx
89                 AtmospCO2(i,j,bi,bj)=dic_pCO2
90               ENDDO
91            ENDDO            ENDDO
92           ENDDO           ENDDO
   
          ENDDO  
93          ENDDO          ENDDO
94        endif        ENDIF
   
 c read from a file (note:  
 c                   gchem_int2=number entries to read  
 c                   gchem_int3=start timestep,  
 c                   gchem_int4=timestep between file entries)  
       if (gchem_int1.eq.2) then  
         if (istate.eq.0) then  
           OPEN(28,FILE='co2atmos.dat',STATUS='old')  
           do it=1,gchem_int2  
             READ(28,*) co2atmos(it)  
             print*,'co2atmos',co2atmos(it)  
           enddo  
          endif  
 c linearly interpolate between file entries  
           ntim=int((myIter-gchem_int3)/gchem_int4)+1  
           aWght = FLOAT(myIter-gchem_int3)  
           bWght = FLOAT(gchem_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  
           tmp=co2atmos(ntim)*bWght+co2atmos(ntim+1)*aWght  
 c         print*,'weights',ntim, aWght, bWght, tmp  
   
           DO bj=myByLo(myThid),myByHi(myThid)  
            DO bi=myBxLo(myThid),myBxHi(myThid)  
95    
96              DO j=1-OLy,sNy+OLy  C read from a file (note:
97    C                   dic_int2=number entries to read
98    C                   dic_int3=start timestep,
99    C                   dic_int4=timestep between file entries)
100          IF (dic_int1.EQ.2) THEN
101    C linearly interpolate between file entries
102            ntim=int((myIter-dic_int3)/dic_int4)+1
103            aWght = FLOAT(myIter-dic_int3)
104            bWght = FLOAT(dic_int4)
105            aWght = 0.5 _d 0 + aWght/bWght - FLOAT(ntim-1)
106            IF (aWght.GT.1. _d 0) THEN
107              ntim=ntim+1
108              aWght=aWght-1. _d 0
109            ENDIF
110            bWght = 1. _d 0 - aWght
111            tmp=co2atmos(ntim)*bWght+co2atmos(ntim+1)*aWght
112            WRITE(ioUnit,*) 'weights',ntim, aWght, bWght, tmp
113    
114            DO bj=myByLo(myThid),myByHi(myThid)
115             DO bi=myBxLo(myThid),myBxHi(myThid)
116               DO j=1-OLy,sNy+OLy
117               DO i=1-OLx,sNx+OLx               DO i=1-OLx,sNx+OLx
   
118                 AtmospCO2(i,j,bi,bj)=tmp                 AtmospCO2(i,j,bi,bj)=tmp
119               ENDDO               ENDDO
             ENDDO  
   
            print*,'AtmospCO2(20,20)',AtmospCO2(20,20,bi,bj)  
   
120             ENDDO             ENDDO
121            ENDDO           ENDDO
122            ENDDO
   
       endif  
   
123    
124  c interactive atmosphere        ENDIF
       if (gchem_int1.eq.3) then  
125    
126  c      _BEGIN_MASTER(myThid)  C interactive atmosphere
127          IF (dic_int1.EQ.3) THEN
128    
129  cMass dry atmosphere = (5.1352+/-0.0003)d18 kg (Trenberth & Smith,  C Mass dry atmosphere = (5.1352+/-0.0003)d18 kg (Trenberth & Smith,
130  cJournal of Climate 2005)  C Journal of Climate 2005)
131  cand Mean molecular mass air = 28.97 g/mol (NASA earth fact sheet)  C and Mean molecular mass air = 28.97 g/mol (NASA earth fact sheet)
132         total_atmos_moles= 1.77 _d 20         total_atmos_moles= 1.77 _d 20
133  c for 278ppmv we need total_atmos_carbon=4.9206e+16  C for 278ppmv we need total_atmos_carbon=4.9206e+16
   
        if (istate.gt.0) then  
         total_ocean_carbon_old=total_ocean_carbon  
         total_atmos_carbon_old=total_atmos_carbon  
        else  
         total_ocean_carbon_old=0. _d 0  
         total_atmos_carbon_old=0. _d 0  
        endif  
   
        total_flux= 0. _d 0  
        total_ocean_carbon= 0. _d 0  
134    
135         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
136          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
137           DO i=1,sNx           tile_flux(bi,bj)   = 0.
138             tile_carbon(bi,bj) = 0.
139             IF (istate.GT.0) THEN
140              DO j=1,sNy
141               DO i=1,sNx
142                 tile_flux(bi,bj) = tile_flux(bi,bj)
143         &                        + FluxCO2(i,j,bi,bj)*rA(i,j,bi,bj)
144         &                         *maskC(i,j,1,bi,bj)*dTtracerLev(1)
145               ENDDO
146              ENDDO
147             ENDIF
148             DO k=1,Nr
149            DO j=1,sNy            DO j=1,sNy
150              if (istate.gt.0) then             DO i=1,sNx
151                total_flux=total_flux+FluxCO2(i,j,bi,bj)*rA(i,j,bi,bj)*               tile_carbon(bi,bj) = tile_carbon(bi,bj)
152       &                        hFacC(i,j,1,bi,bj)*dTtracerLev(1)       &            + ( pTracer(i,j,k,bi,bj,1)
             endif  
             DO k=1,nR  
               total_ocean_carbon= total_ocean_carbon+  
      &              ( Ptracer(i,j,k,bi,bj,1)  
153  #ifdef DIC_BIOTIC  #ifdef DIC_BIOTIC
154       &               +R_cp*Ptracer(i,j,k,bi,bj,4)       &               +R_cp*pTracer(i,j,k,bi,bj,4)
155  #endif  #endif
156       &              ) * rA(i,j,bi,bj)*       &              ) * rA(i,j,bi,bj)
157       &                drF(k)*hFacC(i,j,k,bi,bj)       &                *drF(k)*hFacC(i,j,k,bi,bj)
158              ENDDO             ENDDO
159            ENDDO            ENDDO
160           ENDDO           ENDDO
161          ENDDO          ENDDO
162         ENDDO         ENDDO
163    
164         _GLOBAL_SUM_R8(total_flux,myThid)         CALL GLOBAL_SUM_TILE_RL( tile_flux,   total_flux,   myThid )
165         _GLOBAL_SUM_R8(total_ocean_carbon,myThid)         CALL GLOBAL_SUM_TILE_RL( tile_carbon, total_carbon, myThid )
   
        if (istate.eq.0) then  
 c read state from output file  
           DO i = 1,MAX_LEN_FNAM  
                 fn(i:i) = ' '  
           ENDDO  
           WRITE(fn,'(A,I10.10)') 'dic_atmos.',myIter  
 C     Going to really do some IO. Make everyone except master thread wait.  
            _BARRIER  
 c read in values from last pickup  
            open(26,file=fn,status='old')  
            read(26,*) total_atmos_carbon, atpco2  
            close(26)  
   
        else  
 c calculate new atmos pCO2  
           total_atmos_carbon=total_atmos_carbon - total_flux  
           atpco2=total_atmos_carbon/total_atmos_moles  
 c write out if time for a new pickup  
           permCheckPoint = .FALSE.  
           permCheckPoint =  
      &      DIFFERENT_MULTIPLE(pChkptFreq,myTime,deltaTClock)  
           if (permCheckPoint) then  
             DO i = 1,MAX_LEN_FNAM  
                 fn(i:i) = ' '  
             ENDDO  
             WRITE(fn,'(A,I10.10)') 'dic_atmos.',myIter  
 C     Going to really do some IO. Make everyone except master thread wait.  
             _BARRIER  
 c write values to new pickup  
             open(26,file=fn,status='new')  
             write(26,*) total_atmos_carbon, atpco2  
             close(26)  
166    
167            endif         IF (istate.EQ.0) THEN
168         endif  C use dic_pCO2 as initial atmospheric pCO2 (not restart case):
169             _BEGIN_MASTER(myThid)
170             atpco2 = dic_pCO2
171        atpco2=total_atmos_carbon/total_atmos_moles           total_atmos_carbon = total_atmos_moles*dic_pCO2
172             _END_MASTER(myThid)
173             IF ( nIter0.GT.PTRACERS_Iter0 .OR.
174         &       (nIter0.EQ.PTRACERS_Iter0 .AND. pickupSuff.NE.' ')
175         &      ) THEN
176    C restart case: read previous atmospheric CO2 content & pCO2 from pickup file
177               CALL DIC_READ_CO2_PICKUP( nIter0, myThid )
178             ENDIF
179             _BEGIN_MASTER(myThid)
180    C store initial content:
181             total_ocean_carbon_start=total_carbon
182             total_atmos_carbon_start=total_atmos_carbon
183             total_ocean_carbon_old = total_carbon
184             total_atmos_carbon_old = total_atmos_carbon
185             _END_MASTER(myThid)
186           ELSE
187             _BEGIN_MASTER(myThid)
188    #ifdef ALLOW_AUTODIFF_TAMC
189             atpco2 = dic_pCO2
190    #endif
191    C store previous content:
192             total_ocean_carbon_old = total_ocean_carbon
193             total_atmos_carbon_old = total_atmos_carbon
194    C calculate new atmos pCO2
195             total_atmos_carbon = total_atmos_carbon - total_flux
196             _END_MASTER(myThid)
197           ENDIF
198           _BEGIN_MASTER(myThid)
199           total_ocean_carbon = total_carbon
200           atpco2 = total_atmos_carbon/total_atmos_moles
201    
202           WRITE(ioUnit,*) 'QQ atmos C, total, pCo2',
203         &                     total_atmos_carbon, atpco2
204           total_carbon=total_atmos_carbon + total_ocean_carbon
205           total_carbon_old=total_atmos_carbon_old + total_ocean_carbon_old
206           carbon_diff=total_carbon-total_carbon_old
207           WRITE(ioUnit,*) 'QQ total C, current, old, diff',
208         &                     total_carbon, total_carbon_old, carbon_diff
209           carbon_diff=total_ocean_carbon-total_ocean_carbon_old
210           WRITE(ioUnit,*) 'QQ ocean C, current, old, diff',
211         &         total_ocean_carbon, total_ocean_carbon_old, carbon_diff
212           WRITE(ioUnit,*) 'QQ air-sea flux, addition diff',
213         &                     total_flux, carbon_diff-total_flux
214    
215    C if end of forcing cycle, find total change in ocean carbon
216           IF (istate.EQ.0) THEN
217            total_ocean_carbon_year = total_ocean_carbon
218            total_atmos_carbon_year = total_atmos_carbon
219           ELSE
220            timeCO2budget =
221         &      DIFFERENT_MULTIPLE(externForcingCycle,myTime,deltaTClock)
222            IF ( timeCO2budget ) THEN
223              year_diff_ocean = total_ocean_carbon-total_ocean_carbon_year
224              year_diff_atmos = total_atmos_carbon-total_atmos_carbon_year
225              year_total = (total_ocean_carbon+total_atmos_carbon) -
226         &               (total_ocean_carbon_year+total_atmos_carbon_year)
227              start_diff_ocean = total_ocean_carbon-total_ocean_carbon_start
228              start_diff_atmos = total_atmos_carbon-total_atmos_carbon_start
229              start_total = (total_ocean_carbon+total_atmos_carbon) -
230         &               (total_ocean_carbon_start+total_atmos_carbon_start)
231              WRITE(ioUnit,*) 'QQ YEAR END'
232              WRITE(ioUnit,*) 'year diff: ocean, atmos, total',
233         &                  year_diff_ocean,  year_diff_atmos,  year_total
234              WRITE(ioUnit,*) 'start diff: ocean, atmos, total ',
235         &                 start_diff_ocean, start_diff_atmos, start_total
236    
237              total_ocean_carbon_year = total_ocean_carbon
238              total_atmos_carbon_year = total_atmos_carbon
239            ENDIF
240           ENDIF
241    
242  c     print*,'QQpCO2', total_atmos_carbon, atpco2, total_ocean_carbon,         _END_MASTER(myThid)
243  c    &                 total_flux         _BARRIER
244    
245    C--    Set AtmospCO2 for next iteration:
246         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
247          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
   
248           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
249            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
250               AtmospCO2(i,j,bi,bj)=atpco2              AtmospCO2(i,j,bi,bj) = atpco2
251            ENDDO            ENDDO
252           ENDDO           ENDDO
   
          ENDDO  
253          ENDDO          ENDDO
254           ENDDO
255    
256        print*,'QQ atmos C, total, pCo2', total_atmos_carbon, atpco2        ENDIF
       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  
       print*,'QQ total C, current, old, diff', total_carbon,  
      &                         total_carbon_old, carbon_diff  
       carbon_diff=total_ocean_carbon-total_ocean_carbon_old  
       tmp=carbon_diff-total_flux  
       print*,'QQ ocean C, current, old, diff',total_ocean_carbon,  
      &                   total_ocean_carbon_old, carbon_diff  
       print*,'QQ air-sea flux, addition diff', total_flux, tmp  
   
 c if end of forcing cycle, find total change in ocean carbon  
       if (istate.eq.0) then  
        total_ocean_carbon_start=total_ocean_carbon  
        total_ocean_carbon_year=total_ocean_carbon  
        total_atmos_carbon_start=total_atmos_carbon  
        total_atmos_carbon_year=total_atmos_carbon  
       else  
         permCheckPoint = .FALSE.  
         permCheckPoint =  
      &      DIFFERENT_MULTIPLE(externForcingCycle,myTime,deltaTClock)  
         if (permCheckPoint) 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)  
           print*,'QQ YEAR END'  
           print*,'year diff: ocean, atmos, total', year_diff_ocean,  
      &                year_diff_atmos, year_total  
           print*,'start diff: ocean, atmos, total ', start_diff_ocean,  
      &                start_diff_atmos, start_total  
 c  
           total_ocean_carbon_year=total_ocean_carbon  
           total_atmos_carbon_year=total_atmos_carbon  
         endif  
         endif  
   
 c     _END_MASTER(myThid)  
257    
258        endif  #endif /* ndef USE_ATMOSCO2 */
259    
260  #endif  #endif /* ALLOW_DIC */
 #endif  
261    
262        RETURN        RETURN
263        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22