/[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.3 - (show annotations) (download)
Tue Nov 6 15:54:37 2007 UTC (17 years, 8 months ago) by stephd
Branch: MAIN
Changes since 1.2: +7 -0 lines
o use better estimate atmospheric 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 total_atmos_moles= 1.5 _d 20
158
159 if (istate.gt.0) then
160 total_ocean_carbon_old=total_ocean_carbon
161 total_atmos_carbon_old=total_atmos_carbon
162 else
163 total_ocean_carbon_old=0. _d 0
164 total_atmos_carbon_old=0. _d 0
165 endif
166
167 total_flux= 0. _d 0
168 total_ocean_carbon= 0. _d 0
169
170 DO bj=myByLo(myThid),myByHi(myThid)
171 DO bi=myBxLo(myThid),myBxHi(myThid)
172 DO i=1,sNx
173 DO j=1,sNy
174 if (istate.gt.0) then
175 total_flux=total_flux+FluxCO2(i,j,bi,bj)*rA(i,j,bi,bj)*
176 & hFacC(i,j,1,bi,bj)*dTtracerLev(1)
177 endif
178 DO k=1,nR
179 total_ocean_carbon= total_ocean_carbon+
180 & ( Ptracer(i,j,k,bi,bj,1)+
181 & R_cp*Ptracer(i,j,k,bi,bj,4) )*rA(i,j,bi,bj)*
182 & drF(k)*hFacC(i,j,k,bi,bj)
183 ENDDO
184 ENDDO
185 ENDDO
186 ENDDO
187 ENDDO
188
189 _GLOBAL_SUM_R8(total_flux,myThid)
190 _GLOBAL_SUM_R8(total_ocean_carbon,myThid)
191
192 if (istate.eq.0) then
193 c read state from output file
194 DO i = 1,MAX_LEN_FNAM
195 fn(i:i) = ' '
196 ENDDO
197 WRITE(fn,'(A,I10.10)') 'dic_atmos.',myIter
198 C Going to really do some IO. Make everyone except master thread wait.
199 _BARRIER
200 c read in values from last pickup
201 open(26,file=fn,status='old')
202 read(26,*) total_atmos_carbon, atpco2
203 close(26)
204
205 else
206 c calculate new atmos pCO2
207 total_atmos_carbon=total_atmos_carbon - total_flux
208 atpco2=total_atmos_carbon/total_atmos_moles
209 c write out if time for a new pickup
210 permCheckPoint = .FALSE.
211 permCheckPoint =
212 & DIFFERENT_MULTIPLE(pChkptFreq,myTime,deltaTClock)
213 if (permCheckPoint) then
214 DO i = 1,MAX_LEN_FNAM
215 fn(i:i) = ' '
216 ENDDO
217 WRITE(fn,'(A,I10.10)') 'dic_atmos.',myIter
218 C Going to really do some IO. Make everyone except master thread wait.
219 _BARRIER
220 c write values to new pickup
221 open(26,file=fn,status='new')
222 write(26,*) total_atmos_carbon, atpco2
223 close(26)
224
225 endif
226 endif
227
228
229 atpco2=total_atmos_carbon/total_atmos_moles
230
231 c print*,'QQpCO2', total_atmos_carbon, atpco2, total_ocean_carbon,
232 c & total_flux
233
234 DO bj=myByLo(myThid),myByHi(myThid)
235 DO bi=myBxLo(myThid),myBxHi(myThid)
236
237 DO j=1-OLy,sNy+OLy
238 DO i=1-OLx,sNx+OLx
239 AtmospCO2(i,j,bi,bj)=atpco2
240 ENDDO
241 ENDDO
242
243 ENDDO
244 ENDDO
245
246 print*,'QQ atmos C, total, pCo2', total_atmos_carbon, atpco2
247 total_carbon=total_atmos_carbon + total_ocean_carbon
248 total_carbon_old=total_atmos_carbon_old + total_ocean_carbon_old
249 carbon_diff=total_carbon-total_carbon_old
250 print*,'QQ total C, current, old, diff', total_carbon,
251 & total_carbon_old, carbon_diff
252 carbon_diff=total_ocean_carbon-total_ocean_carbon_old
253 tmp=carbon_diff-total_flux
254 print*,'QQ ocean C, current, old, diff',total_ocean_carbon,
255 & total_ocean_carbon_old, carbon_diff
256 print*,'QQ air-sea flux, addition diff', total_flux, tmp
257
258 c if end of forcing cycle, find total change in ocean carbon
259 if (istate.eq.0) then
260 total_ocean_carbon_start=total_ocean_carbon
261 total_ocean_carbon_year=total_ocean_carbon
262 total_atmos_carbon_start=total_atmos_carbon
263 total_atmos_carbon_year=total_atmos_carbon
264 else
265 permCheckPoint = .FALSE.
266 permCheckPoint =
267 & DIFFERENT_MULTIPLE(externForcingCycle,myTime,deltaTClock)
268 if (permCheckPoint) then
269 year_diff_ocean=total_ocean_carbon-total_ocean_carbon_year
270 year_diff_atmos=total_atmos_carbon-total_atmos_carbon_year
271 year_total=(total_ocean_carbon+total_atmos_carbon) -
272 & (total_ocean_carbon_year+total_atmos_carbon_year)
273 start_diff_ocean=total_ocean_carbon-total_ocean_carbon_start
274 start_diff_atmos=total_atmos_carbon-total_atmos_carbon_start
275 start_total=(total_ocean_carbon+total_atmos_carbon) -
276 & (total_ocean_carbon_start+total_atmos_carbon_start)
277 print*,'QQ YEAR END'
278 print*,'year diff: ocean, atmos, total', year_diff_ocean,
279 & year_diff_atmos, year_total
280 print*,'start diff: ocean, atmos, total ', start_diff_ocean,
281 & start_diff_atmos, start_total
282 c
283 total_ocean_carbon_year=total_ocean_carbon
284 total_atmos_carbon_year=total_atmos_carbon
285 endif
286 endif
287
288 c _END_MASTER(myThid)
289
290 endif
291
292 #endif
293 #endif
294
295 RETURN
296 END

  ViewVC Help
Powered by ViewVC 1.1.22