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 -- |