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

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

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


Revision 1.5 - (hide annotations) (download)
Thu Nov 8 22:35:23 2007 UTC (17 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59j
Changes since 1.4: +23 -13 lines
- add cvs $Header:$ and $Name:$
- no time to fix multi-threaded now: put a stop.
- change argument list order: myThid in last position
- standard way to compute coeff. for time interpolation (aWght is computed
  directly in real*8 instead of in real*4 and then converted)

1 jmc 1.5 C $Header: $
2     C $Name: $
3    
4 stephd 1.1 #include "DIC_OPTIONS.h"
5     #include "PTRACERS_OPTIONS.h"
6     #include "GCHEM_OPTIONS.h"
7    
8     CBOP
9     C !ROUTINE: DIC_ATMOS
10    
11     C !INTERFACE: ==========================================================
12 jmc 1.5 SUBROUTINE DIC_ATMOS( istate, myTime, myIter, myThid )
13 stephd 1.1
14     C !DESCRIPTION:
15     C Calculate the atmospheric pCO2
16     C gchem_int1:
17     C 0=use default 278.d-6
18 jmc 1.5 C 1=use constant value - gchem_rl1, read in from data.gchem
19     C 2=read in from file
20 stephd 1.1 C 3=interact with atmospheric box
21     C !USES: ===============================================================
22     IMPLICIT NONE
23     #include "SIZE.h"
24     #include "DYNVARS.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27     #include "GRID.h"
28     #include "FFIELDS.h"
29     #include "DIC_ABIOTIC.h"
30     #ifdef DIC_BIOTIC
31     #include "PTRACERS_SIZE.h"
32     #include "PTRACERS_FIELDS.h"
33     #include "DIC_BIOTIC.h"
34     #endif
35     #include "GCHEM.h"
36     #include "DIC_ATMOS.h"
37    
38     C !INPUT PARAMETERS: ===================================================
39     C myThid :: thread number
40     C myIter :: current timestep
41     C myTime :: current time
42     C istate :: 0=initial call, 1=subsequent calls
43     INTEGER myIter, myThid, istate
44     _RL myTime
45    
46     #ifdef ALLOW_PTRACERS
47     LOGICAL DIFFERENT_MULTIPLE
48     EXTERNAL DIFFERENT_MULTIPLE
49    
50     C !LOCAL VARIABLES: ====================================================
51     INTEGER bi, bj, I,J,k
52     INTEGER it, ntim
53     c
54     _RL total_flux
55     _RL total_ocean_carbon_old
56     _RL total_atmos_carbon_old
57     _RL total_atmos_moles
58     _RL atpco2
59     _RL total_carbon_old, total_carbon, carbon_diff
60     _RL tmp
61     _RL year_diff_ocean, year_diff_atmos, year_total
62     _RL start_diff_ocean, start_diff_atmos, start_total
63     C variables for reading CO2 input files
64     _RL aWght, bWght
65     c
66     CHARACTER*(MAX_LEN_FNAM) fn
67     LOGICAL permCheckPoint
68     CEOP
69    
70     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
71    
72     c if coupled to atmsopheric model, use the
73     c Co2 value passed from the coupler
74     #ifndef USE_ATMOSCO2
75    
76 jmc 1.5 IF ( nThreads .GT. 1 ) THEN
77     C Problem with I/O and global-sum for multi-threaded execution
78     C Needs to be fixed before using this S/R in multi-threaded run
79     STOP 'S/R DIC_ATMOS: multi-threaded not right'
80     ENDIF
81    
82 stephd 1.1 c default - set only once
83     if (gchem_int1.eq.0.and.istate.eq.0) then
84     DO bj=myByLo(myThid),myByHi(myThid)
85     DO bi=myBxLo(myThid),myBxHi(myThid)
86    
87     DO j=1-OLy,sNy+OLy
88     DO i=1-OLx,sNx+OLx
89     AtmospCO2(i,j,bi,bj)=278.0 _d -6
90     ENDDO
91     ENDDO
92    
93     ENDDO
94     ENDDO
95     endif
96    
97     c user specified value - set only once
98     if (gchem_int1.eq.1.and.istate.eq.0) then
99     DO bj=myByLo(myThid),myByHi(myThid)
100     DO bi=myBxLo(myThid),myBxHi(myThid)
101    
102     DO j=1-OLy,sNy+OLy
103     DO i=1-OLx,sNx+OLx
104     AtmospCO2(i,j,bi,bj)=gchem_rl1
105     ENDDO
106     ENDDO
107    
108     ENDDO
109     ENDDO
110     endif
111    
112     c read from a file (note:
113     c gchem_int2=number entries to read
114     c gchem_int3=start timestep,
115     c gchem_int4=timestep between file entries)
116     if (gchem_int1.eq.2) then
117 jmc 1.5 if (istate.eq.0) then
118 stephd 1.1 OPEN(28,FILE='co2atmos.dat',STATUS='old')
119     do it=1,gchem_int2
120     READ(28,*) co2atmos(it)
121     print*,'co2atmos',co2atmos(it)
122     enddo
123     endif
124     c linearly interpolate between file entries
125     ntim=int((myIter-gchem_int3)/gchem_int4)+1
126 jmc 1.5 aWght = FLOAT(myIter-gchem_int3)
127     bWght = FLOAT(gchem_int4)
128     aWght = 0.5 _d 0 + aWght/bWght - FLOAT(ntim-1)
129     if (aWght.gt.1. _d 0) then
130 stephd 1.1 ntim=ntim+1
131 jmc 1.5 aWght=aWght-1. _d 0
132 stephd 1.1 endif
133 jmc 1.5 bWght = 1. _d 0 - aWght
134 stephd 1.1 tmp=co2atmos(ntim)*bWght+co2atmos(ntim+1)*aWght
135     c print*,'weights',ntim, aWght, bWght, tmp
136    
137     DO bj=myByLo(myThid),myByHi(myThid)
138     DO bi=myBxLo(myThid),myBxHi(myThid)
139    
140     DO j=1-OLy,sNy+OLy
141     DO i=1-OLx,sNx+OLx
142    
143     AtmospCO2(i,j,bi,bj)=tmp
144     ENDDO
145     ENDDO
146    
147     print*,'AtmospCO2(20,20)',AtmospCO2(20,20,bi,bj)
148 jmc 1.5
149 stephd 1.1 ENDDO
150     ENDDO
151    
152    
153     endif
154    
155    
156     c interactive atmosphere
157     if (gchem_int1.eq.3) then
158    
159 stephd 1.2 c _BEGIN_MASTER(myThid)
160 stephd 1.1
161 stephd 1.3 cMass dry atmosphere = (5.1352+/-0.0003)d18 kg (Trenberth & Smith,
162     cJournal of Climate 2005)
163     cand Mean molecular mass air = 28.97 g/mol (NASA earth fact sheet)
164     total_atmos_moles= 1.77 _d 20
165     c for 278ppmv we need total_atmos_carbon=4.9206e+16
166    
167 stephd 1.1 if (istate.gt.0) then
168     total_ocean_carbon_old=total_ocean_carbon
169     total_atmos_carbon_old=total_atmos_carbon
170     else
171     total_ocean_carbon_old=0. _d 0
172     total_atmos_carbon_old=0. _d 0
173     endif
174    
175     total_flux= 0. _d 0
176     total_ocean_carbon= 0. _d 0
177    
178     DO bj=myByLo(myThid),myByHi(myThid)
179     DO bi=myBxLo(myThid),myBxHi(myThid)
180     DO i=1,sNx
181     DO j=1,sNy
182     if (istate.gt.0) then
183     total_flux=total_flux+FluxCO2(i,j,bi,bj)*rA(i,j,bi,bj)*
184     & hFacC(i,j,1,bi,bj)*dTtracerLev(1)
185 jmc 1.5 endif
186 stephd 1.1 DO k=1,nR
187     total_ocean_carbon= total_ocean_carbon+
188     & ( Ptracer(i,j,k,bi,bj,1)+
189     & R_cp*Ptracer(i,j,k,bi,bj,4) )*rA(i,j,bi,bj)*
190     & drF(k)*hFacC(i,j,k,bi,bj)
191     ENDDO
192     ENDDO
193     ENDDO
194     ENDDO
195     ENDDO
196 stephd 1.2
197 stephd 1.1 _GLOBAL_SUM_R8(total_flux,myThid)
198     _GLOBAL_SUM_R8(total_ocean_carbon,myThid)
199    
200     if (istate.eq.0) then
201     c read state from output file
202     DO i = 1,MAX_LEN_FNAM
203     fn(i:i) = ' '
204     ENDDO
205     WRITE(fn,'(A,I10.10)') 'dic_atmos.',myIter
206     C Going to really do some IO. Make everyone except master thread wait.
207     _BARRIER
208     c read in values from last pickup
209     open(26,file=fn,status='old')
210     read(26,*) total_atmos_carbon, atpco2
211     close(26)
212    
213     else
214     c calculate new atmos pCO2
215     total_atmos_carbon=total_atmos_carbon - total_flux
216     atpco2=total_atmos_carbon/total_atmos_moles
217     c write out if time for a new pickup
218     permCheckPoint = .FALSE.
219     permCheckPoint =
220     & DIFFERENT_MULTIPLE(pChkptFreq,myTime,deltaTClock)
221     if (permCheckPoint) then
222     DO i = 1,MAX_LEN_FNAM
223     fn(i:i) = ' '
224     ENDDO
225     WRITE(fn,'(A,I10.10)') 'dic_atmos.',myIter
226     C Going to really do some IO. Make everyone except master thread wait.
227     _BARRIER
228 stephd 1.3 c write values to new pickup
229 stephd 1.1 open(26,file=fn,status='new')
230     write(26,*) total_atmos_carbon, atpco2
231     close(26)
232    
233     endif
234     endif
235    
236    
237     atpco2=total_atmos_carbon/total_atmos_moles
238    
239 jmc 1.5 c print*,'QQpCO2', total_atmos_carbon, atpco2, total_ocean_carbon,
240 stephd 1.1 c & total_flux
241    
242     DO bj=myByLo(myThid),myByHi(myThid)
243     DO bi=myBxLo(myThid),myBxHi(myThid)
244    
245     DO j=1-OLy,sNy+OLy
246     DO i=1-OLx,sNx+OLx
247     AtmospCO2(i,j,bi,bj)=atpco2
248     ENDDO
249     ENDDO
250    
251     ENDDO
252     ENDDO
253    
254     print*,'QQ atmos C, total, pCo2', total_atmos_carbon, atpco2
255     total_carbon=total_atmos_carbon + total_ocean_carbon
256     total_carbon_old=total_atmos_carbon_old + total_ocean_carbon_old
257     carbon_diff=total_carbon-total_carbon_old
258 jmc 1.5 print*,'QQ total C, current, old, diff', total_carbon,
259 stephd 1.1 & total_carbon_old, carbon_diff
260     carbon_diff=total_ocean_carbon-total_ocean_carbon_old
261     tmp=carbon_diff-total_flux
262     print*,'QQ ocean C, current, old, diff',total_ocean_carbon,
263     & total_ocean_carbon_old, carbon_diff
264     print*,'QQ air-sea flux, addition diff', total_flux, tmp
265    
266     c if end of forcing cycle, find total change in ocean carbon
267     if (istate.eq.0) then
268     total_ocean_carbon_start=total_ocean_carbon
269     total_ocean_carbon_year=total_ocean_carbon
270     total_atmos_carbon_start=total_atmos_carbon
271     total_atmos_carbon_year=total_atmos_carbon
272     else
273     permCheckPoint = .FALSE.
274     permCheckPoint =
275     & DIFFERENT_MULTIPLE(externForcingCycle,myTime,deltaTClock)
276     if (permCheckPoint) then
277     year_diff_ocean=total_ocean_carbon-total_ocean_carbon_year
278     year_diff_atmos=total_atmos_carbon-total_atmos_carbon_year
279     year_total=(total_ocean_carbon+total_atmos_carbon) -
280     & (total_ocean_carbon_year+total_atmos_carbon_year)
281     start_diff_ocean=total_ocean_carbon-total_ocean_carbon_start
282     start_diff_atmos=total_atmos_carbon-total_atmos_carbon_start
283     start_total=(total_ocean_carbon+total_atmos_carbon) -
284     & (total_ocean_carbon_start+total_atmos_carbon_start)
285     print*,'QQ YEAR END'
286     print*,'year diff: ocean, atmos, total', year_diff_ocean,
287     & year_diff_atmos, year_total
288     print*,'start diff: ocean, atmos, total ', start_diff_ocean,
289     & start_diff_atmos, start_total
290     c
291     total_ocean_carbon_year=total_ocean_carbon
292     total_atmos_carbon_year=total_atmos_carbon
293     endif
294     endif
295    
296 stephd 1.2 c _END_MASTER(myThid)
297 stephd 1.1
298     endif
299    
300     #endif
301     #endif
302    
303     RETURN
304     END

  ViewVC Help
Powered by ViewVC 1.1.22