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

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

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


Revision 1.4 - (show annotations) (download)
Tue Nov 6 15:57:45 2007 UTC (17 years, 8 months ago) by stephd
Branch: MAIN
Changes since 1.3: +0 -2 lines
o remove old estimate of atmoshperic moles!!

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

  ViewVC Help
Powered by ViewVC 1.1.22