/[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.1 - (show annotations) (download)
Mon Nov 5 21:01:37 2007 UTC (17 years, 8 months ago) by stephd
Branch: MAIN
o adding flexibility for determining atmospheric pCO2
  (now defaults to 278d-6, can be read in from data.gchem, can be
   read in from file for time-varing values, and can be linked
   to interactive well mixed atmospheric box)

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

  ViewVC Help
Powered by ViewVC 1.1.22