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