/[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.5 - (show 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 C $Header: $
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_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 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 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 if (istate.eq.0) then
118 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 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 ntim=ntim+1
131 aWght=aWght-1. _d 0
132 endif
133 bWght = 1. _d 0 - aWght
134 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
149 ENDDO
150 ENDDO
151
152
153 endif
154
155
156 c interactive atmosphere
157 if (gchem_int1.eq.3) then
158
159 c _BEGIN_MASTER(myThid)
160
161 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 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 endif
186 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
197 _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 c write values to new pickup
229 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 c print*,'QQpCO2', total_atmos_carbon, atpco2, total_ocean_carbon,
240 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 print*,'QQ total C, current, old, diff', total_carbon,
259 & 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 c _END_MASTER(myThid)
297
298 endif
299
300 #endif
301 #endif
302
303 RETURN
304 END

  ViewVC Help
Powered by ViewVC 1.1.22