/[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.8 - (show annotations) (download)
Fri Apr 4 21:32:37 2008 UTC (17 years, 3 months ago) by dfer
Branch: MAIN
Changes since 1.7: +2 -5 lines
Merging DIC_ABIOTIC.h and DIC_BIOTIC.h

1 C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_atmos.F,v 1.7 2008/04/04 18:57:36 dfer Exp $
2 C $Name: $
3
4 #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 SUBROUTINE DIC_ATMOS( istate, myTime, myIter, myThid )
13
14 C !DESCRIPTION:
15 C Calculate the atmospheric pCO2
16 C gchem_int1:
17 C 0=use default 278.d-6
18 C 1=use constant value - gchem_rl1, read in from data.gchem
19 C 2=read in from file
20 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_VARS.h"
30 #include "PTRACERS_SIZE.h"
31 #include "PTRACERS_FIELDS.h"
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 IF ( nThreads .GT. 1 .AND.
74 & ( gchem_int1.EQ.2 .OR. gchem_int1.EQ.3 ) ) THEN
75 C Problem with I/O and global-sum for multi-threaded execution
76 C Needs to be fixed before using this S/R in multi-threaded run
77 STOP 'S/R DIC_ATMOS: multi-threaded not right'
78 ENDIF
79
80 c default - set only once
81 if (gchem_int1.eq.0.and.istate.eq.0) then
82 DO bj=myByLo(myThid),myByHi(myThid)
83 DO bi=myBxLo(myThid),myBxHi(myThid)
84
85 DO j=1-OLy,sNy+OLy
86 DO i=1-OLx,sNx+OLx
87 AtmospCO2(i,j,bi,bj)=278.0 _d -6
88 ENDDO
89 ENDDO
90
91 ENDDO
92 ENDDO
93 endif
94
95 c user specified value - set only once
96 if (gchem_int1.eq.1.and.istate.eq.0) then
97 DO bj=myByLo(myThid),myByHi(myThid)
98 DO bi=myBxLo(myThid),myBxHi(myThid)
99
100 DO j=1-OLy,sNy+OLy
101 DO i=1-OLx,sNx+OLx
102 AtmospCO2(i,j,bi,bj)=gchem_rl1
103 ENDDO
104 ENDDO
105
106 ENDDO
107 ENDDO
108 endif
109
110 c read from a file (note:
111 c gchem_int2=number entries to read
112 c gchem_int3=start timestep,
113 c gchem_int4=timestep between file entries)
114 if (gchem_int1.eq.2) then
115 if (istate.eq.0) then
116 OPEN(28,FILE='co2atmos.dat',STATUS='old')
117 do it=1,gchem_int2
118 READ(28,*) co2atmos(it)
119 print*,'co2atmos',co2atmos(it)
120 enddo
121 endif
122 c linearly interpolate between file entries
123 ntim=int((myIter-gchem_int3)/gchem_int4)+1
124 aWght = FLOAT(myIter-gchem_int3)
125 bWght = FLOAT(gchem_int4)
126 aWght = 0.5 _d 0 + aWght/bWght - FLOAT(ntim-1)
127 if (aWght.gt.1. _d 0) then
128 ntim=ntim+1
129 aWght=aWght-1. _d 0
130 endif
131 bWght = 1. _d 0 - aWght
132 tmp=co2atmos(ntim)*bWght+co2atmos(ntim+1)*aWght
133 c print*,'weights',ntim, aWght, bWght, tmp
134
135 DO bj=myByLo(myThid),myByHi(myThid)
136 DO bi=myBxLo(myThid),myBxHi(myThid)
137
138 DO j=1-OLy,sNy+OLy
139 DO i=1-OLx,sNx+OLx
140
141 AtmospCO2(i,j,bi,bj)=tmp
142 ENDDO
143 ENDDO
144
145 print*,'AtmospCO2(20,20)',AtmospCO2(20,20,bi,bj)
146
147 ENDDO
148 ENDDO
149
150
151 endif
152
153
154 c interactive atmosphere
155 if (gchem_int1.eq.3) then
156
157 c _BEGIN_MASTER(myThid)
158
159 cMass dry atmosphere = (5.1352+/-0.0003)d18 kg (Trenberth & Smith,
160 cJournal of Climate 2005)
161 cand Mean molecular mass air = 28.97 g/mol (NASA earth fact sheet)
162 total_atmos_moles= 1.77 _d 20
163 c for 278ppmv we need total_atmos_carbon=4.9206e+16
164
165 if (istate.gt.0) then
166 total_ocean_carbon_old=total_ocean_carbon
167 total_atmos_carbon_old=total_atmos_carbon
168 else
169 total_ocean_carbon_old=0. _d 0
170 total_atmos_carbon_old=0. _d 0
171 endif
172
173 total_flux= 0. _d 0
174 total_ocean_carbon= 0. _d 0
175
176 DO bj=myByLo(myThid),myByHi(myThid)
177 DO bi=myBxLo(myThid),myBxHi(myThid)
178 DO i=1,sNx
179 DO j=1,sNy
180 if (istate.gt.0) then
181 total_flux=total_flux+FluxCO2(i,j,bi,bj)*rA(i,j,bi,bj)*
182 & hFacC(i,j,1,bi,bj)*dTtracerLev(1)
183 endif
184 DO k=1,nR
185 total_ocean_carbon= total_ocean_carbon+
186 & ( Ptracer(i,j,k,bi,bj,1)
187 #ifdef DIC_BIOTIC
188 & +R_cp*Ptracer(i,j,k,bi,bj,4)
189 #endif
190 & ) * rA(i,j,bi,bj)*
191 & drF(k)*hFacC(i,j,k,bi,bj)
192 ENDDO
193 ENDDO
194 ENDDO
195 ENDDO
196 ENDDO
197
198 _GLOBAL_SUM_R8(total_flux,myThid)
199 _GLOBAL_SUM_R8(total_ocean_carbon,myThid)
200
201 if (istate.eq.0) then
202 c read state from output file
203 DO i = 1,MAX_LEN_FNAM
204 fn(i:i) = ' '
205 ENDDO
206 WRITE(fn,'(A,I10.10)') 'dic_atmos.',myIter
207 C Going to really do some IO. Make everyone except master thread wait.
208 _BARRIER
209 c read in values from last pickup
210 open(26,file=fn,status='old')
211 read(26,*) total_atmos_carbon, atpco2
212 close(26)
213
214 else
215 c calculate new atmos pCO2
216 total_atmos_carbon=total_atmos_carbon - total_flux
217 atpco2=total_atmos_carbon/total_atmos_moles
218 c write out if time for a new pickup
219 permCheckPoint = .FALSE.
220 permCheckPoint =
221 & DIFFERENT_MULTIPLE(pChkptFreq,myTime,deltaTClock)
222 if (permCheckPoint) then
223 DO i = 1,MAX_LEN_FNAM
224 fn(i:i) = ' '
225 ENDDO
226 WRITE(fn,'(A,I10.10)') 'dic_atmos.',myIter
227 C Going to really do some IO. Make everyone except master thread wait.
228 _BARRIER
229 c write values to new pickup
230 open(26,file=fn,status='new')
231 write(26,*) total_atmos_carbon, atpco2
232 close(26)
233
234 endif
235 endif
236
237
238 atpco2=total_atmos_carbon/total_atmos_moles
239
240 c print*,'QQpCO2', total_atmos_carbon, atpco2, total_ocean_carbon,
241 c & total_flux
242
243 DO bj=myByLo(myThid),myByHi(myThid)
244 DO bi=myBxLo(myThid),myBxHi(myThid)
245
246 DO j=1-OLy,sNy+OLy
247 DO i=1-OLx,sNx+OLx
248 AtmospCO2(i,j,bi,bj)=atpco2
249 ENDDO
250 ENDDO
251
252 ENDDO
253 ENDDO
254
255 print*,'QQ atmos C, total, pCo2', total_atmos_carbon, atpco2
256 total_carbon=total_atmos_carbon + total_ocean_carbon
257 total_carbon_old=total_atmos_carbon_old + total_ocean_carbon_old
258 carbon_diff=total_carbon-total_carbon_old
259 print*,'QQ total C, current, old, diff', total_carbon,
260 & total_carbon_old, carbon_diff
261 carbon_diff=total_ocean_carbon-total_ocean_carbon_old
262 tmp=carbon_diff-total_flux
263 print*,'QQ ocean C, current, old, diff',total_ocean_carbon,
264 & total_ocean_carbon_old, carbon_diff
265 print*,'QQ air-sea flux, addition diff', total_flux, tmp
266
267 c if end of forcing cycle, find total change in ocean carbon
268 if (istate.eq.0) then
269 total_ocean_carbon_start=total_ocean_carbon
270 total_ocean_carbon_year=total_ocean_carbon
271 total_atmos_carbon_start=total_atmos_carbon
272 total_atmos_carbon_year=total_atmos_carbon
273 else
274 permCheckPoint = .FALSE.
275 permCheckPoint =
276 & DIFFERENT_MULTIPLE(externForcingCycle,myTime,deltaTClock)
277 if (permCheckPoint) then
278 year_diff_ocean=total_ocean_carbon-total_ocean_carbon_year
279 year_diff_atmos=total_atmos_carbon-total_atmos_carbon_year
280 year_total=(total_ocean_carbon+total_atmos_carbon) -
281 & (total_ocean_carbon_year+total_atmos_carbon_year)
282 start_diff_ocean=total_ocean_carbon-total_ocean_carbon_start
283 start_diff_atmos=total_atmos_carbon-total_atmos_carbon_start
284 start_total=(total_ocean_carbon+total_atmos_carbon) -
285 & (total_ocean_carbon_start+total_atmos_carbon_start)
286 print*,'QQ YEAR END'
287 print*,'year diff: ocean, atmos, total', year_diff_ocean,
288 & year_diff_atmos, year_total
289 print*,'start diff: ocean, atmos, total ', start_diff_ocean,
290 & start_diff_atmos, start_total
291 c
292 total_ocean_carbon_year=total_ocean_carbon
293 total_atmos_carbon_year=total_atmos_carbon
294 endif
295 endif
296
297 c _END_MASTER(myThid)
298
299 endif
300
301 #endif
302 #endif
303
304 RETURN
305 END

  ViewVC Help
Powered by ViewVC 1.1.22