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

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

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


Revision 1.1 - (show annotations) (download)
Wed Jun 25 21:00:36 2003 UTC (20 years, 11 months ago) by stephd
Branch: MAIN
CVS Tags: checkpoint51a_post, checkpoint51c_post, checkpoint51b_pre, checkpoint51b_post
initial checking in biogeochemistry packages

1 cswdcost -- add sunroutine ---
2 #include "CPP_OPTIONS.h"
3 #include "GCHEM_OPTIONS.h"
4
5
6 CStartOfInterFace
7 SUBROUTINE DIC_DIAGS(
8 I myTime,myIter,myThid)
9
10 C /==========================================================\
11 C | SUBROUTINE DIC_COST i |
12 C |==========================================================|
13 IMPLICIT NONE
14
15 C == GLobal variables ==
16 #include "SIZE.h"
17 #include "DYNVARS.h"
18 #include "EEPARAMS.h"
19 #include "PARAMS.h"
20 #include "GRID.h"
21 #include "PTRACERS.h"
22 #include "GCHEM.h"
23 #include "DIC_ABIOTIC.h"
24 #ifdef DIC_BIOTIC
25 #include "DIC_BIOTIC.h"
26 #include "DIC_DIAGS.h"
27 #endif
28 #ifdef ALLOW_SEAICE
29 #include "ICE.h"
30 #endif
31
32 C == Routine arguments ==
33 INTEGER myIter
34 _RL myTime
35 INTEGER myThid
36
37 #ifdef ALLOW_DIC_COST
38
39 C == Local variables ==
40 LOGICAL DIFFERENT_MULTIPLE
41 EXTERNAL DIFFERENT_MULTIPLE
42 INTEGER i, j, bi, bj, k, it
43 _RL po4av(nR)
44 _RL o2av(nR)
45 _RL volvar(nR)
46 cswdmonth -add-
47 _RL po4avm(12,4)
48 _RL o2avm(12,4)
49 _RL rdt
50 INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
51
52 cswddmonth -- end-
53 c
54 c initialize to zero
55 IF ( myIter .EQ. nIter0+1 ) THEN
56 DO bj = myByLo(myThid), myByHi(myThid)
57 DO bi = myBxLo(myThid), myBxHi(myThid)
58 CALL TIMEAVE_RESET(PO4ann, Nr, bi, bj, myThid)
59 CALL TIMEAVE_RESET(O2ann, Nr, bi, bj, myThid)
60 CALL TIMEAVE_RESET(PO4obs, Nr, bi, bj, myThid)
61 CALL TIMEAVE_RESET(O2obs, Nr, bi, bj, myThid)
62 cswdmonth
63 CALL TIMEAVE_RESET(PO4obsl1, Nr, bi, bj, myThid)
64 CALL TIMEAVE_RESET(PO4obsl2, Nr, bi, bj, myThid)
65 CALL TIMEAVE_RESET(PO4obsl3, Nr, bi, bj, myThid)
66 cQQ CALL TIMEAVE_RESET(PO4obsl4, Nr, bi, bj, myThid)
67 CALL TIMEAVE_RESET(O2obsl1, Nr, bi, bj, myThid)
68 CALL TIMEAVE_RESET(O2obsl2, Nr, bi, bj, myThid)
69 CALL TIMEAVE_RESET(O2obsl3, Nr, bi, bj, myThid)
70 cQQ CALL TIMEAVE_RESET(O2obsl4, Nr, bi, bj, myThid)
71 CALL TIMEAVE_RESET(PO4lev1, Nr, bi, bj, myThid)
72 CALL TIMEAVE_RESET(PO4lev2, Nr, bi, bj, myThid)
73 CALL TIMEAVE_RESET(PO4lev3, Nr, bi, bj, myThid)
74 cQQ CALL TIMEAVE_RESET(PO4lev4, Nr, bi, bj, myThid)
75 CALL TIMEAVE_RESET(O2lev1, Nr, bi, bj, myThid)
76 CALL TIMEAVE_RESET(O2lev2, Nr, bi, bj, myThid)
77 cQQ CALL TIMEAVE_RESET(O2lev3, Nr, bi, bj, myThid)
78 cQQ CALL TIMEAVE_RESET(O2lev4, Nr, bi, bj, myThid)
79 cswdmonth -end-
80 do k=1,Nr
81 OBS_Timetave(bi,bj,k)=0.d0
82 po4av(k)=0.d0
83 o2av(k)=0.d0
84 po4var(k)=0.d0
85 o2var(k)=0.d0
86 volvar(k)=0.d0
87 enddo
88 cswdmonth
89 do k=1,3
90 do it=1,12
91 OBSM_Timetave(bi,bj,it)=0.d0
92 po4avm(it,k)=0.d0
93 o2avm(it,k)=0.d0
94 po4varm(it,k)=0.d0
95 o2varm(it,k)=0.d0
96 enddo
97 enddo
98 ENDDO
99 ENDDO
100 _BEGIN_MASTER( myThid )
101 CALL READ_FLD_XYZ_RL( 'input/po4obs.bin', ' ',
102 & po4obs, 0, myThid )
103 CALL READ_FLD_XYZ_RL( 'input/o2obs.bin', ' ',
104 & o2obs, 0, myThid )
105 cswdmonth
106 CALL READ_FLD_XYZ_RL( 'input/po4lev1.bin', ' ',
107 & po4obsl1, 0, myThid )
108 CALL READ_FLD_XYZ_RL( 'input/po4lev2.bin', ' ',
109 & po4obsl2, 0, myThid )
110 CALL READ_FLD_XYZ_RL( 'input/po4lev3.bin', ' ',
111 & po4obsl3, 0, myThid )
112 cQQ CALL READ_FLD_XYZ_RL( 'input/po4lev4.bin', ' ',
113 cQQ & po4obsl4, 0, myThid )
114 CALL READ_FLD_XYZ_RL( 'input/o2lev1.bin', ' ',
115 & o2obsl1, 0, myThid )
116 CALL READ_FLD_XYZ_RL( 'input/o2lev2.bin', ' ',
117 & o2obsl2, 0, myThid )
118 CALL READ_FLD_XYZ_RL( 'input/o2lev3.bin', ' ',
119 & o2obsl3, 0, myThid )
120 cQQ CALL READ_FLD_XYZ_RL( 'input/o2lev4.bin', ' ',
121 cQQ & o2obsl4, 0, myThid )
122 cswdmonth -end-
123 _END_MASTER(myThid)
124 _EXCH_XYZ_R8(po4obs , myThid )
125 _EXCH_XYZ_R8(o2obs , myThid )
126 cswdmonth -add-
127 _EXCH_XYZ_R8(po4obsl1 , myThid )
128 _EXCH_XYZ_R8(po4obsl2 , myThid )
129 _EXCH_XYZ_R8(po4obsl3 , myThid )
130 cQQ _EXCH_XYZ_R8(po4obsl4 , myThid )
131 _EXCH_XYZ_R8(o2obsl1 , myThid )
132 _EXCH_XYZ_R8(o2obsl2 , myThid )
133 _EXCH_XYZ_R8(o2obsl3 , myThid )
134 cQQ _EXCH_XYZ_R8(o2obsl4 , myThid )
135 cswdmonth -end-
136
137 _BARRIER
138 c calculate layer means
139 _BEGIN_MASTER( mythid )
140 do k=1,Nr
141 call tracer_meanarea(myIter,myTime,myThid,po4obs, k,
142 & po4av(k))
143 call tracer_meanarea(myIter,myTime,myThid,o2obs, k,
144 & o2av(k))
145 c print*,po4av(k), o2av(k)
146 enddo
147 cswdmonth -add-
148 do it=1,12
149 call tracer_meanarea(myIter,myTime,myThid,po4obsl1,it,
150 & po4avm(it,1))
151 call tracer_meanarea(myIter,myTime,myThid,po4obsl2,it,
152 & po4avm(it,2))
153 call tracer_meanarea(myIter,myTime,myThid,po4obsl3,it,
154 & po4avm(it,3))
155 cQQ call tracer_meanarea(myIter,myTime,myThid,po4obsl4,it,
156 cQQ & po4avm(it,4))
157 call tracer_meanarea(myIter,myTime,myThid,o2obsl1,it,
158 & o2avm(it,1))
159 call tracer_meanarea(myIter,myTime,myThid,o2obsl2,it,
160 & o2avm(it,2))
161 call tracer_meanarea(myIter,myTime,myThid,o2obsl3,it,
162 & o2avm(it,3))
163 cQQ call tracer_meanarea(myIter,myTime,myThid,o2obsl4,it,
164 cQQ & o2avm(it,4))
165
166 enddo
167 _END_MASTER(myThid)
168 c calculate layer variance
169 _BEGIN_MASTER( mythid )
170 DO bj = myByLo(myThid), myByHi(myThid)
171 DO bi = myBxLo(myThid), myBxHi(myThid)
172 DO j=1-OLy,sNy+OLy
173 DO i=1-OLx,sNx+OLx
174 DO k=1,Nr
175 volvar(k)=volvar(k)+
176 & rA(i,j,bi,bj)*drF(k)*maskC(i,j,k,bi,bj)
177 po4var(k)=po4var(k)+(po4obs(i,j,k,bi,bj)-po4av(k))**2
178 & *rA(i,j,bi,bj)*drF(k)*maskC(i,j,k,bi,bj)
179 o2var(k)=o2var(k)+(o2obs(i,j,k,bi,bj)-o2av(k))**2
180 & *rA(i,j,bi,bj)*drF(k)*maskC(i,j,k,bi,bj)
181 ENDDO
182 cswdmonth -add-
183 DO it=1,12
184 po4varm(it,1)=po4varm(it,1)+
185 & (po4obsl1(i,j,it,bi,bj)-po4avm(it,1))**2
186 & *rA(i,j,bi,bj)*drF(1)*maskC(i,j,1,bi,bj)
187 po4varm(it,2)=po4varm(it,2)+
188 & (po4obsl2(i,j,it,bi,bj)-po4avm(it,2))**2
189 & *rA(i,j,bi,bj)*drF(2)*maskC(i,j,2,bi,bj)
190 po4varm(it,3)=po4varm(it,3)+
191 & (po4obsl3(i,j,it,bi,bj)-po4avm(it,3))**2
192 & *rA(i,j,bi,bj)*drF(3)*maskC(i,j,3,bi,bj)
193 cQQ po4varm(it,4)=po4varm(it,4)+
194 cQQ & (po4obsl4(i,j,it,bi,bj)-po4avm(it,4))**2
195 cQQ & *rA(i,j,bi,bj)*drF(4)*maskC(i,j,4,bi,bj)
196 o2varm(it,1)=o2varm(it,1)+
197 & (o2obsl1(i,j,it,bi,bj)-o2avm(it,1))**2
198 & *rA(i,j,bi,bj)*drF(1)*maskC(i,j,1,bi,bj)
199 o2varm(it,2)=o2varm(it,2)+
200 & (o2obsl2(i,j,it,bi,bj)-o2avm(it,2))**2
201 & *rA(i,j,bi,bj)*drF(2)*maskC(i,j,2,bi,bj)
202 o2varm(it,3)=o2varm(it,3)+
203 & (o2obsl3(i,j,it,bi,bj)-o2avm(it,3))**2
204 & *rA(i,j,bi,bj)*drF(3)*maskC(i,j,3,bi,bj)
205 cQQ o2varm(it,4)=o2varm(it,4)+
206 cQQ & (o2obsl4(i,j,it,bi,bj)-o2avm(it,4))**2
207 cQQ & *rA(i,j,bi,bj)*drF(4)*maskC(i,j,4,bi,bj)
208
209 ENDDO
210 ENDDO
211 ENDDO
212 ENDDO
213 ENDDO
214 DO k=1,Nr
215 po4var(k)=po4var(k)/volvar(k)
216 o2var(k)=o2var(k)/volvar(k)
217 cQQ print*,po4var(k),o2var(k)
218 ENDDO
219 cswdmonth- add-
220 DO k=1,3
221 Do it=1,12
222 po4varm(it,k)=po4varm(it,k)/volvar(k)
223 o2varm(it,k)=o2varm(it,k)/volvar(k)
224 ENDDO
225 ENDDO
226 cswdmonth -end-
227 _END_MASTER(myThid)
228 ENDIF
229 C
230 c averages
231 DO bj = myByLo(myThid), myByHi(myThid)
232 DO bi = myBxLo(myThid), myBxHi(myThid)
233 DO k=1,Nr
234 OBS_timetave(bi,bj,k)=OBS_timetave(bi,bj,k)+
235 & deltaTclock
236 DO j=1-OLy,sNy+OLy
237 DO i=1-OLx,sNx+OLx
238 po4ann(i,j,k,bi,bj)=po4ann(i,j,k,bi,bj)+
239 & PTRACER(i,j,k,bi,bj,3)*deltaTclock
240 o2ann(i,j,k,bi,bj)=o2ann(i,j,k,bi,bj)+
241 & PTRACER(i,j,k,bi,bj,5)*deltaTclock
242 ENDDO
243 ENDDO
244 ENDDO
245 ENDDO
246 ENDDO
247 cswdmonth-add--
248 rdt=1. _d 0 / deltaTclock
249 nForcingPeriods=int(externForcingCycle/externForcingPeriod+0.5)
250 Imytm=int(myTime*rdt+0.5)
251 Ifprd=int(externForcingPeriod*rdt+0.5)
252 Ifcyc=int(externForcingCycle*rdt+0.5)
253 Iftm=mod( Imytm+Ifcyc ,Ifcyc)
254 it=int(Iftm/Ifprd)+1
255 c print*,'QQ timing check', mytime, myIter, it
256 DO bj = myByLo(myThid), myByHi(myThid)
257 DO bi = myBxLo(myThid), myBxHi(myThid)
258 OBSM_timetave(bi,bj,it)=OBSM_timetave(bi,bj,it)+
259 & deltaTclock
260 DO j=1-OLy,sNy+OLy
261 DO i=1-OLx,sNx+OLx
262 po4lev1(i,j,it,bi,bj)=po4lev1(i,j,it,bi,bj)+
263 & PTRACER(i,j,1,bi,bj,3)*deltaTclock
264 po4lev2(i,j,it,bi,bj)=po4lev2(i,j,it,bi,bj)+
265 & PTRACER(i,j,2,bi,bj,3)*deltaTclock
266 po4lev3(i,j,it,bi,bj)=po4lev3(i,j,it,bi,bj)+
267 & PTRACER(i,j,3,bi,bj,3)*deltaTclock
268 cQQ po4lev4(i,j,it,bi,bj)=po4lev4(i,j,it,bi,bj)+
269 cQQ & PTRACER(i,j,4,bi,bj,3)*deltaTclock
270 o2lev1(i,j,it,bi,bj)=o2lev1(i,j,it,bi,bj)+
271 & PTRACER(i,j,1,bi,bj,5)*deltaTclock
272 o2lev2(i,j,it,bi,bj)=o2lev2(i,j,it,bi,bj)+
273 & PTRACER(i,j,2,bi,bj,5)*deltaTclock
274 o2lev3(i,j,it,bi,bj)=o2lev3(i,j,it,bi,bj)+
275 & PTRACER(i,j,3,bi,bj,5)*deltaTclock
276 cQQ O2lev4(i,j,it,bi,bj)=O2lev4(i,j,it,bi,bj)+
277 cQQ & PTRACER(i,j,4,bi,bj,5)*deltaTclock
278 ENDDO
279 ENDDO
280 ENDDO
281 ENDDO
282
283 cswdmonth-end--
284
285
286 C Dump files and restart average computation if needed
287 IF ( DIFFERENT_MULTIPLE(
288 & externForcingCycle,myTime,myTime-deltaTClock)
289 & ) THEN
290 C
291 C Normalize by integrated time
292 DO bj = myByLo(myThid), myByHi(myThid)
293 DO bi = myBxLo(myThid), myBxHi(myThid)
294 CALL TIMEAVE_NORMALIZ(PO4ann, OBS_timetave, nR ,
295 & bi,bj,myThid)
296 CALL TIMEAVE_NORMALIZ(O2ann, OBS_timetave, nR ,
297 & bi,bj,myThid)
298 cswdmonth-add
299 CALL TIMEAVE_NORMALIZ(PO4lev1, OBSM_timetave, 12 ,
300 & bi,bj,myThid)
301 CALL TIMEAVE_NORMALIZ(PO4lev2, OBSM_timetave, 12 ,
302 & bi,bj,myThid)
303 CALL TIMEAVE_NORMALIZ(PO4lev3, OBSM_timetave, 12 ,
304 & bi,bj,myThid)
305 cQQ CALL TIMEAVE_NORMALIZ(PO4lev4, OBSM_timetave, 12 ,
306 cQQ & bi,bj,myThid)
307 CALL TIMEAVE_NORMALIZ(O2lev1, OBSM_timetave, 12 ,
308 & bi,bj,myThid)
309 CALL TIMEAVE_NORMALIZ(O2lev2, OBSM_timetave, 12 ,
310 & bi,bj,myThid)
311 CALL TIMEAVE_NORMALIZ(O2lev3, OBSM_timetave, 12 ,
312 & bi,bj,myThid)
313 cQQ CALL TIMEAVE_NORMALIZ(O2lev4, OBSM_timetave, 12 ,
314 cQQ & bi,bj,myThid)
315 cswdmonth -end-
316 ENDDO
317 ENDDO
318 c
319 c calculate cost function
320 _BEGIN_MASTER( mythid )
321 call DIC_COST(myTime,myIter,myThid)
322 _END_MASTER( mythid )
323 c
324 C Reset averages to zero
325 DO bj = myByLo(myThid), myByHi(myThid)
326 DO bi = myBxLo(myThid), myBxHi(myThid)
327 CALL TIMEAVE_RESET(PO4ann,Nr,bi,bj,myThid)
328 CALL TIMEAVE_RESET(O2ann,Nr,bi,bj,myThid)
329 do k=1,Nr
330 OBS_Timetave(bi,bj,k)=0.d0
331 enddo
332 cswdmonth-add--
333 CALL TIMEAVE_RESET(PO4lev1,12,bi,bj,myThid)
334 CALL TIMEAVE_RESET(PO4lev2,12,bi,bj,myThid)
335 CALL TIMEAVE_RESET(PO4lev3,12,bi,bj,myThid)
336 cQQ CALL TIMEAVE_RESET(PO4lev4,12,bi,bj,myThid)
337 CALL TIMEAVE_RESET(O2lev1,12,bi,bj,myThid)
338 CALL TIMEAVE_RESET(O2lev2,12,bi,bj,myThid)
339 CALL TIMEAVE_RESET(O2lev3,12,bi,bj,myThid)
340 cQQ CALL TIMEAVE_RESET(O2lev4,12,bi,bj,myThid)
341 do it=1,12
342 OBSM_Timetave(bi,bj,it)=0.d0
343 enddo
344 cswdmonth-end--
345 ENDDO
346 ENDDO
347 c
348 ENDIF
349 #endif
350 c
351 RETURN
352 END
353 cswd -- end added subroutine --

  ViewVC Help
Powered by ViewVC 1.1.22