/[MITgcm]/MITgcm_contrib/arctic40km/code_ad/ecco_cost_weights.F
ViewVC logotype

Annotation of /MITgcm_contrib/arctic40km/code_ad/ecco_cost_weights.F

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

Revision 1.2 - (hide annotations) (download)
Wed May 31 21:30:36 2006 UTC (19 years, 1 month ago) by heimbach
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +307 -199 lines 
o adding input_ad
o various changes to adjoint setup compared to forward,
  some of which for initial testing
1 heimbach 1.2 C $Header: /u/gcmpack/MITgcm/pkg/ecco/ecco_cost_weights.F,v 1.17 2006/05/12 13:40:18 heimbach Exp $
2 heimbach 1.1
3     #include "COST_CPPOPTIONS.h"
4    
5    
6     subroutine ecco_cost_weights( mythid )
7    
8     c ==================================================================
9     c SUBROUTINE ecco_cost_weights
10     c ==================================================================
11     c
12     c o Read the weights used for the cost function evaluation.
13     c
14     c started: Christian Eckert eckert@mit.edu 30-Jun-1999
15     c
16     c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
17     c
18     c - Restructured the code in order to create a package
19     c for the MITgcmUV.
20     c
21     c Christian Eckert eckert@mit.edu 02-May-2000
22     c
23     c - corrected typo in mdsreadfield( sflux_errfile );
24     c wp --> wsflux. Spotted by Patrick Heimbach.
25     c
26     c ==================================================================
27     c SUBROUTINE ecco_cost_weights
28     c ==================================================================
29    
30     implicit none
31    
32     c == global variables ==
33    
34     #include "EEPARAMS.h"
35     #include "SIZE.h"
36     #include "PARAMS.h"
37     #include "GRID.h"
38    
39     #include "ctrl.h"
40     #include "ecco_cost.h"
41    
42     c == routine arguments ==
43    
44     integer mythid
45    
46     c == local variables ==
47    
48     integer bi,bj
49     integer i,j,k
50     integer itlo,ithi
51     integer jtlo,jthi
52     integer jmin,jmax
53     integer imin,imax
54     integer gwunit
55     integer irec,nnz
56     integer ilo,ihi
57     integer iobcs
58    
59     _RL factor
60     _RL wti(nr)
61     _RL wsi(nr)
62     _RL wui(nr)
63     _RL wvi(nr)
64     _RL whflux0m
65     _RL wsflux0m
66     _RL wtau0m
67     _RL ratio
68     _RL dummy
69    
70     c == external ==
71    
72     integer ifnblnk
73     external ifnblnk
74     integer ilnblnk
75     external ilnblnk
76    
77     c == end of interface ==
78    
79     jtlo = mybylo(mythid)
80     jthi = mybyhi(mythid)
81     itlo = mybxlo(mythid)
82     ithi = mybxhi(mythid)
83     jmin = 1-oly
84     jmax = sny+oly
85     imin = 1-olx
86     imax = snx+olx
87    
88     c-- Initialize background weights
89 heimbach 1.2 whflux0m = whflux0
90     wsflux0m = wsflux0
91     wtau0m = wtau0
92 heimbach 1.1
93     c-- Initialize variance (weight) fields.
94     do k = 1,nr
95     wti(k) = 0. _d 0
96     wsi(k) = 0. _d 0
97     wui(k) = 0. _d 0
98     wvi(k) = 0. _d 0
99     enddo
100     do bj = jtlo,jthi
101     do bi = itlo,ithi
102     do j = jmin,jmax
103     do i = imin,imax
104     whflux (i,j,bi,bj) = 0. _d 0
105     whfluxm (i,j,bi,bj) = 0. _d 0
106     wsflux (i,j,bi,bj) = 0. _d 0
107     wsfluxm (i,j,bi,bj) = 0. _d 0
108     wtauu (i,j,bi,bj) = 0. _d 0
109     wtauum (i,j,bi,bj) = 0. _d 0
110     wtauv (i,j,bi,bj) = 0. _d 0
111     wtauvm (i,j,bi,bj) = 0. _d 0
112     watemp (i,j,bi,bj) = 0. _d 0
113     waqh (i,j,bi,bj) = 0. _d 0
114 heimbach 1.2 wprecip (i,j,bi,bj) = 0. _d 0
115     wswflux (i,j,bi,bj) = 0. _d 0
116     wswdown (i,j,bi,bj) = 0. _d 0
117 heimbach 1.1 wuwind (i,j,bi,bj) = 0. _d 0
118     wvwind (i,j,bi,bj) = 0. _d 0
119     wsst (i,j,bi,bj) = 0. _d 0
120     wsss (i,j,bi,bj) = 0. _d 0
121     wtp (i,j,bi,bj) = 0. _d 0
122     wers (i,j,bi,bj) = 0. _d 0
123     wp (i,j,bi,bj) = 0. _d 0
124     wudrift (i,j,bi,bj) = 0. _d 0
125     wvdrift (i,j,bi,bj) = 0. _d 0
126     cph(
127 heimbach 1.2 whflux2 (i,j,bi,bj) = 0. _d 0
128     wsflux2 (i,j,bi,bj) = 0. _d 0
129     wtauu2 (i,j,bi,bj) = 0. _d 0
130     wtauv2 (i,j,bi,bj) = 0. _d 0
131 heimbach 1.1 cph)
132     enddo
133     enddo
134     enddo
135     enddo
136     do bj = jtlo,jthi
137     do bi = itlo,ithi
138     do k = 1,Nr
139     wtheta (k,bi,bj) = 0. _d 0
140     wsalt (k,bi,bj) = 0. _d 0
141     wctdt (k,bi,bj) = 0. _d 0
142     wctds (k,bi,bj) = 0. _d 0
143 heimbach 1.2 wdiffkr(k,bi,bj) = wdiffkr0
144     wkapgm (k,bi,bj) = wkapgm0
145     wedtaux(k,bi,bj) = wedtau0
146     wedtauy(k,bi,bj) = wedtau0
147 heimbach 1.1 do j = jmin,jmax
148     do i = imin,imax
149     wtheta2 (i,j,k,bi,bj) = 0. _d 0
150     wsalt2 (i,j,k,bi,bj) = 0. _d 0
151 heimbach 1.2 wdiffkr2(i,j,k,bi,bj) = wdiffkr0
152     wkapgm2 (i,j,k,bi,bj) = wkapgm0
153     wedtaux2(i,j,k,bi,bj) = wedtau0
154     wedtauy2(i,j,k,bi,bj) = wedtau0
155 heimbach 1.1 wthetaLev (i,j,k,bi,bj) = 0. _d 0
156     wsaltLev (i,j,k,bi,bj) = 0. _d 0
157 heimbach 1.2 wdiffkrFld(i,j,k,bi,bj) = wdiffkr0
158     wkapgmFld (i,j,k,bi,bj) = wkapgm0
159     wedtauxFld(i,j,k,bi,bj) = wedtau0
160     wedtauyFld(i,j,k,bi,bj) = wedtau0
161 heimbach 1.1 enddo
162     enddo
163     enddo
164     enddo
165     enddo
166    
167     #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || \
168     defined (ALLOW_OBCS_CONTROL))
169     do iobcs = 1,nobcs
170     do k = 1,Nr
171     #if (defined (ALLOW_OBCSN_CONTROL) || \
172     defined (ALLOW_OBCSN_COST_CONTRIBUTION))
173     wobcsn(k,iobcs) = 0. _d 0
174     #endif
175     #if (defined (ALLOW_OBCSS_CONTROL) || \
176     defined (ALLOW_OBCSS_COST_CONTRIBUTION))
177     wobcss(k,iobcs) = 0. _d 0
178     #endif
179     #if (defined (ALLOW_OBCSW_CONTROL) || \
180     defined (ALLOW_OBCSW_COST_CONTRIBUTION))
181     wobcsw(k,iobcs) = 0. _d 0
182     #endif
183     #if (defined (ALLOW_OBCSE_CONTROL) || \
184     defined (ALLOW_OBCSE_COST_CONTRIBUTION))
185     wobcse(k,iobcs) = 0. _d 0
186     #endif
187     enddo
188     enddo
189     #endif
190    
191     c-- Build area weighting matrix used in the cost function
192     c-- contributions.
193    
194     c-- Define frame.
195     do j = jmin,jmax
196     do i = imin,imax
197     c-- North/South and West/East edges set to zero.
198     if ( (j .lt. 1) .or. (j .gt. sny) .or.
199     & (i .lt. 1) .or. (i .gt. snx) ) then
200     frame(i,j) = 0. _d 0
201     else
202     frame(i,j) = 1. _d 0
203     endif
204     enddo
205     enddo
206    
207     c-- First account for the grid used.
208     if (usingCartesianGrid) then
209     factor = 0. _d 0
210     else if (usingSphericalPolarGrid) then
211     factor = 1. _d 0
212     endif
213    
214     do bj = jtlo,jthi
215     do bi = itlo,ithi
216     do j = jmin,jmax
217     do i = imin,imax
218     cds cosphi(i,j,bi,bj) = cos(yc(i,j,bi,bj)*deg2rad*factor)*
219     cds & frame(i,j)
220     cosphi(i,j,bi,bj) = frame(i,j)
221     enddo
222     enddo
223     enddo
224     enddo
225    
226     c-- Read error information and set up weight matrices.
227     _BEGIN_MASTER(myThid)
228     ilo = ifnblnk(data_errfile)
229     ihi = ilnblnk(data_errfile)
230 heimbach 1.2 CALL OPEN_COPY_DATA_FILE(
231 heimbach 1.1 I data_errfile(ilo:ihi),
232     I 'ECCO_COST_WEIGHTS',
233     O gwunit,
234     I myThid )
235    
236 heimbach 1.2 read(gwunit,*) ratio
237 heimbach 1.1 #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || defined (ALLOW_OBCS_CONTROL))
238     & , wbaro
239     #endif
240     do k = 1,nr
241     read(gwunit,*) wti(k), wsi(k)
242 heimbach 1.2 #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || defined (ALLOW_OBCS_CONTROL))
243 heimbach 1.1 & , wvi(k)
244     #endif
245     end do
246     close(gwunit)
247    
248     _END_MASTER(myThid)
249    
250     _BARRIER
251    
252     cph jmin = 1
253     cph jmax = sny
254     cph imin = 1
255     cph imax = snx
256    
257     do bj = jtlo,jthi
258     do bi = itlo,ithi
259    
260     wsfluxmm(bi,bj) = 1.
261     whfluxmm(bi,bj) = 1.
262    
263     c-- The "classic" state estimation tool wastes memory here;
264     c-- as long as there is not more information available there
265     c-- is no need to add the zonal and meridional directions.
266     do k = 1,nr
267     wtheta(k,bi,bj) = wti(k)
268     wsalt (k,bi,bj) = wsi(k)
269     wcurrent(k,bi,bj) = wvi(k)
270     enddo
271    
272     do k = 1,nr
273     #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
274     wobcsn(k,1) = wti(k)
275     wobcsn(k,2) = wsi(k)
276 heimbach 1.2 wobcsn(k,3) = wvi(k)
277     wobcsn(k,4) = wvi(k)
278 heimbach 1.1 #endif
279     #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
280     wobcss(k,1) = wti(k)
281     wobcss(k,2) = wsi(k)
282 heimbach 1.2 wobcss(k,3) = wvi(k)
283     wobcss(k,4) = wvi(k)
284 heimbach 1.1 #endif
285     #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
286     wobcsw(k,1) = wti(k)
287     wobcsw(k,2) = wsi(k)
288 heimbach 1.2 wobcsw(k,3) = wvi(k)
289     wobcsw(k,4) = wvi(k)
290 heimbach 1.1 #endif
291     #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
292     wobcse(k,1) = wti(k)
293     wobcse(k,2) = wsi(k)
294 heimbach 1.2 wobcse(k,3) = wvi(k)
295     wobcse(k,4) = wvi(k)
296 heimbach 1.1 #endif
297     enddo
298     enddo
299     enddo
300    
301     #ifdef ALLOW_SALT0_COST_CONTRIBUTION
302     if ( salt0errfile .NE. ' ' ) then
303     call mdsreadfield( salt0errfile, 32, 'RL', Nr,
304     & wsaltLev, 1, mythid)
305     do bj = jtlo,jthi
306     do bi = itlo,ithi
307     do k = 1,nr
308     do j = jmin,jmax
309     do i = imin,imax
310     c-- Test for missing values.
311     if ( wsaltLev(i,j,k,bi,bj).eq.0 ) then
312     wsaltLev(i,j,k,bi,bj) = 0. _d 0
313     else
314 heimbach 1.2 wsaltLev(i,j,k,bi,bj)=frame(i,j)/(
315     $ wsaltLev(i,j,k,bi,bj)*wsaltLev(i,j,k,bi,bj) )
316 heimbach 1.1 endif
317     enddo
318     enddo
319     enddo
320     enddo
321     enddo
322     else
323     do bj = jtlo,jthi
324     do bi = itlo,ithi
325     do k = 1,nr
326     do j = jmin,jmax
327     do i = imin,imax
328     wsaltLev(i,j,k,bi,bj)=ratio/(wsalt(k,bi,bj)
329     $ *wsalt(k,bi,bj)) *frame(i,j)
330     enddo
331     enddo
332     enddo
333     enddo
334 heimbach 1.2 enddo
335 heimbach 1.1 endif
336 heimbach 1.2 call active_write_xyz( 'wsaltLev', wsaltLev,
337 heimbach 1.1 & 1, 0, mythid, dummy)
338     #endif
339    
340     #ifdef ALLOW_THETA0_COST_CONTRIBUTION
341     if ( temp0errfile .NE. ' ' ) then
342     call mdsreadfield( temp0errfile, 32, 'RL', Nr,
343     & wthetaLev, 1, mythid)
344     do bj = jtlo,jthi
345     do bi = itlo,ithi
346     do k = 1,nr
347     do j = jmin,jmax
348     do i = imin,imax
349     c-- Test for missing values.
350     if ( wthetaLev(i,j,k,bi,bj).eq.0 ) then
351     wthetaLev(i,j,k,bi,bj) = 0. _d 0
352     else
353 heimbach 1.2 wthetaLev(i,j,k,bi,bj)=frame(i,j)/(
354     $ wthetaLev(i,j,k,bi,bj)*wthetaLev(i,j,k,bi,bj) )
355 heimbach 1.1 endif
356     enddo
357     enddo
358     enddo
359     enddo
360     enddo
361     else
362     do bj = jtlo,jthi
363     do bi = itlo,ithi
364     do k = 1,nr
365     do j = jmin,jmax
366     do i = imin,imax
367     wthetaLev(i,j,k,bi,bj)=ratio/(wtheta(k,bi,bj)
368     $ *wtheta(k,bi,bj)) *frame(i,j)
369     enddo
370     enddo
371     enddo
372     enddo
373 heimbach 1.2 enddo
374 heimbach 1.1 endif
375 heimbach 1.2 call active_write_xyz( 'wthetaLev', wthetaLev,
376 heimbach 1.1 & 1, 0, mythid, dummy)
377     #endif
378    
379     #if (defined (ALLOW_ARGO_SALT_COST_CONTRIBUTION) || \
380     defined (ALLOW_CTDS_COST_CONTRIBUTION)|| \
381     defined (ALLOW_CTDSCLIM_COST_CONTRIBUTION))
382    
383     if ( salterrfile .NE. ' ' ) then
384     call mdsreadfield( salterrfile, 32, 'RL', Nr, wsalt2, 1,
385     & mythid)
386    
387     do bj = jtlo,jthi
388     do bi = itlo,ithi
389     do k = 1,nr
390     do j = jmin,jmax
391     do i = imin,imax
392     c-- Test for missing values.
393     if (wsalt(k,bi,bj).eq.0 .and.
394     $ wsalt2(i,j,k,bi,bj).eq.0) then
395     wsalt2(i,j,k,bi,bj) = 0. _d 0
396     else
397 heimbach 1.2 wsalt2(i,j,k,bi,bj)=frame(i,j)/
398     $ ( wsalt2(i,j,k,bi,bj)*wsalt2(i,j,k,bi,bj) )
399 heimbach 1.1 endif
400     enddo
401     enddo
402     enddo
403     enddo
404     enddo
405     else
406     do bj = jtlo,jthi
407     do bi = itlo,ithi
408     do k = 1,nr
409     do j = jmin,jmax
410     do i = imin,imax
411     wsalt2(i,j,k,bi,bj)=ratio/(wsalt(k,bi,bj)
412     $ *wsalt(k,bi,bj)) *frame(i,j)
413     enddo
414     enddo
415     enddo
416     enddo
417     enddo
418     endif
419     #endif
420    
421     #if (defined (ALLOW_ARGO_THETA_COST_CONTRIBUTION) || \
422     defined (ALLOW_CTDT_COST_CONTRIBUTION) || \
423     defined (ALLOW_CTDTCLIM_COST_CONTRIBUTION) || \
424     defined (ALLOW_XBT_COST_CONTRIBUTION))
425    
426     if ( temperrfile .NE. ' ' ) then
427     call mdsreadfield( temperrfile, 32, 'RL', Nr, wtheta2, 1,
428     & mythid)
429     do bj = jtlo,jthi
430     do bi = itlo,ithi
431     do k = 1,nr
432     do j = jmin,jmax
433     do i = imin,imax
434     c-- Test for missing values.
435     if (wtheta(k,bi,bj).eq.0 .and.
436     $ wtheta2(i,j,k,bi,bj).eq.0) then
437     wtheta2(i,j,k,bi,bj) = 0. _d 0
438     else
439 heimbach 1.2 wtheta2(i,j,k,bi,bj)= frame(i,j)/
440     $ ( wtheta2(i,j,k,bi,bj)*wtheta2(i,j,k,bi,bj) )
441 heimbach 1.1 endif
442     enddo
443     enddo
444     enddo
445     enddo
446     enddo
447     else
448     do bj = jtlo,jthi
449     do bi = itlo,ithi
450     do k = 1,nr
451     do j = jmin,jmax
452     do i = imin,imax
453     if (wtheta(k,bi,bj).eq.0 ) then
454     wtheta2(i,j,k,bi,bj) = 0. _d 0
455     else
456     wtheta2(i,j,k,bi,bj) = ratio/(wtheta(k,bi,bj)
457     $ *wtheta(k,bi,bj))*frame(i,j)
458     endif
459     enddo
460     enddo
461     enddo
462     enddo
463     enddo
464     endif
465     #endif
466    
467     k = 1
468     do bj = jtlo,jthi
469     do bi = itlo,ithi
470     do j = jmin,jmax
471     do i = imin,imax
472     if (_hFacC(i,j,k,bi,bj) .eq. 0.) then
473     wsst(i,j,bi,bj) = 0. _d 0
474     wsss(i,j,bi,bj) = 0. _d 0
475     else
476     cph wsst(i,j,bi,bj) = wtheta(k,bi,bj)*10.
477     cph wsss(i,j,bi,bj) = wsalt(k,bi,bj)*10.
478     cph factor 5. by D Stammer
479     wsst(i,j,bi,bj) = wtheta(k,bi,bj)*frame(i,j)
480     wsss(i,j,bi,bj) = wsalt(k,bi,bj)*frame(i,j)
481     endif
482     enddo
483     enddo
484     enddo
485     enddo
486     #ifdef ALLOW_EGM96_ERROR_DIAG
487     c-- Read egm-96 geoid covariance. Data in units of meters.
488     nnz = 1
489     irec = 1
490     call mdsreadfield( geoid_errfile, cost_iprec, cost_yftype, nnz,
491     & wp, irec, mythid )
492     c-- Set all tile edges to zero.
493     do bj = jtlo,jthi
494     do bi = itlo,ithi
495     do j = jmin,jmax
496     do i = imin,imax
497     wp(i,j,bi,bj) = wp(i,j,bi,bj)*frame(i,j)
498     cph-indonesian(
499 heimbach 1.2 if ( xC(i,j,bi,bj) .GT. 120. .AND.
500 heimbach 1.1 & xC(i,j,bi,bj) .LT. 130. .AND.
501     & yC(i,j,bi,bj) .GT. -10. .AND.
502     & yC(i,j,bi,bj) .LT. 10. ) then
503 heimbach 1.2 wp(i,j,bi,bj) = wp(i,j,bi,bj)*100.
504 heimbach 1.1 endif
505     cph-indonesian)
506     enddo
507     enddo
508     enddo
509     enddo
510     #else
511     do bj = jtlo,jthi
512     do bi = itlo,ithi
513     do j = jmin,jmax
514     do i = imin,imax
515     wp(i,j,bi,bj) = frame(i,j)
516     enddo
517     enddo
518     enddo
519     enddo
520     #endif
521    
522     #ifdef ALLOW_SSH_COST_CONTRIBUTION
523     c-- Read T/P SSH anomaly rms field. Data in units of centimeters.
524     nnz = 1
525     irec = 1
526     call mdsreadfield( ssh_errfile, cost_iprec, cost_yftype, nnz,
527     & wtp, irec, mythid )
528    
529     do bj = jtlo,jthi
530     do bi = itlo,ithi
531     k = 1
532     do j = jmin,jmax
533     do i = imin,imax
534     c-- Unit conversion to meters. ERS error is set to
535     c-- T/P error + 5cm
536     if (_hFacC(i,j,k,bi,bj) .eq. 0.) then
537     wtp (i,j,bi,bj) = 0. _d 0
538     wers(i,j,bi,bj) = 0. _d 0
539     else
540     wtp (i,j,bi,bj) = ( wtp(i,j,bi,bj) * 0.01 * 0.5 )
541     & *frame(i,j)
542     wers(i,j,bi,bj) = ( wtp(i,j,bi,bj) + 0.05 )
543     & *frame(i,j)
544     endif
545     cph-indonesian(
546 heimbach 1.2 if ( xC(i,j,bi,bj) .GT. 120. .AND.
547 heimbach 1.1 & xC(i,j,bi,bj) .LT. 130. .AND.
548     & yC(i,j,bi,bj) .GT. -10. .AND.
549     & yC(i,j,bi,bj) .LT. 10. ) then
550 heimbach 1.2 wtp(i,j,bi,bj) = wtp(i,j,bi,bj)*100.
551     wers(i,j,bi,bj) = wers(i,j,bi,bj)*100.
552 heimbach 1.1 endif
553     cph-indonesian)
554     enddo
555     enddo
556     enddo
557     enddo
558     #endif /* ALLOW_SSH_COST_CONTRIBUTION */
559    
560     c-- Read zonal wind stress variance.
561     #if (defined (ALLOW_SCAT_COST_CONTRIBUTION))
562    
563     nnz = 1
564     irec = 1
565     call mdsreadfield( scatx_errfile, cost_iprec, cost_yftype, nnz,
566     & wscatx, irec, mythid )
567     call mdsreadfield( scaty_errfile, cost_iprec, cost_yftype, nnz,
568     & wscaty, irec, mythid )
569    
570     do bj = jtlo,jthi
571     do bi = itlo,ithi
572     k = 1
573     do j = jmin,jmax
574     do i = imin,imax
575     c-- Test for missing values.
576     if (wscatx(i,j,bi,bj) .lt. -9900.) then
577     wscatx(i,j,bi,bj) = 0. _d 0
578     endif
579     wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)
580     wscatx(i,j,bi,bj) = max(wscatx(i,j,bi,bj),wtau0)
581 heimbach 1.2 wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)*maskw(i,j,k,bi,bj)
582 heimbach 1.1 & *frame(i,j)
583     if (wscaty(i,j,bi,bj) .lt. -9900.) then
584     wscaty(i,j,bi,bj) = 0. _d 0
585     endif
586     wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)
587     wscaty(i,j,bi,bj) = max(wscaty(i,j,bi,bj),wtau0)
588 heimbach 1.2 wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)*masks(i,j,k,bi,bj)
589 heimbach 1.1 & *frame(i,j)
590     enddo
591     enddo
592     enddo
593     enddo
594    
595     #endif
596    
597     c-- Read zonal wind stress variance.
598     #if (defined (ALLOW_STRESS_MEAN_COST_CONTRIBUTION))
599     nnz = 1
600     irec = 1
601     cph call mdsreadfield( tauum_errfile, cost_iprec, cost_yftype,
602     cph & nnz, wtauum, irec, mythid )
603     cph call mdsreadfield( tauvm_errfile, cost_iprec, cost_yftype,
604     cph & nnz, wtauvm, irec, mythid )
605    
606     do bj = jtlo,jthi
607     do bi = itlo,ithi
608     k = 1
609     do j = jmin,jmax
610     do i = imin,imax
611     c-- Test for missing values.
612     if (wtauum(i,j,bi,bj) .lt. -9900.) then
613     wtauum(i,j,bi,bj) = 0. _d 0
614     endif
615     wtauum(i,j,bi,bj) = max(wtauum(i,j,bi,bj),wtau0m)
616     & *frame(i,j)
617     c-- Test for missing values.
618     if (wtauvm(i,j,bi,bj) .lt. -9900.) then
619     wtauvm(i,j,bi,bj) = 0. _d 0
620     endif
621     wtauvm(i,j,bi,bj) = max(wtauvm(i,j,bi,bj),wtau0m)
622     & *frame(i,j)
623     enddo
624     enddo
625     enddo
626     enddo
627     #endif
628    
629 heimbach 1.2 #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION) || defined (ALLOW_USTRESS_CONTROL))
630 heimbach 1.1 nnz = 1
631     ce irec = 2
632     ce( due to Patrick's processing:
633     irec = 1
634     ce)
635     if ( tauu_errfile .NE. ' ' ) then
636     call mdsreadfield( tauu_errfile, cost_iprec, cost_yftype, nnz,
637 heimbach 1.2 & wtauu, irec, mythid )
638 heimbach 1.1 endif
639    
640     do bj = jtlo,jthi
641     do bi = itlo,ithi
642     k = 1
643     do j = jmin,jmax
644     do i = imin,imax
645     c-- Test for missing values.
646     if (wtauu(i,j,bi,bj) .lt. -9900.) then
647     wtauu(i,j,bi,bj) = 0. _d 0
648     endif
649     wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*0.1
650     wtauu(i,j,bi,bj) = max(wtauu(i,j,bi,bj),wtau0)
651 heimbach 1.2 wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*maskw(i,j,k,bi,bj)
652 heimbach 1.1 & *frame(i,j)
653     cph(
654 heimbach 1.2 wtauu2(i,j,bi,bj) = wtau0*maskW(i,j,k,bi,bj)*frame(i,j)
655 heimbach 1.1 cph)
656     enddo
657     enddo
658     enddo
659     enddo
660    
661 heimbach 1.2 #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION) || defined (ALLOW_UWIND_CONTROL))
662 heimbach 1.1
663     nnz = 1
664     ce irec = 2
665     ce( due to Patrick's processing:
666     irec = 1
667     ce)
668     if ( uwind_errfile .NE. ' ' ) then
669     call mdsreadfield( uwind_errfile, cost_iprec, cost_yftype, nnz,
670 heimbach 1.2 & wuwind, irec, mythid )
671 heimbach 1.1 endif
672    
673     do bj = jtlo,jthi
674     do bi = itlo,ithi
675     k = 1
676     do j = jmin,jmax
677     do i = imin,imax
678     c-- Test for missing values.
679     if (wuwind(i,j,bi,bj) .lt. -9900.) then
680     wuwind(i,j,bi,bj) = 0. _d 0
681     endif
682     wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)
683     wuwind(i,j,bi,bj) = max(wuwind(i,j,bi,bj),wwind0)
684 heimbach 1.2 wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)*maskc(i,j,k,bi,bj)
685 heimbach 1.1 & *frame(i,j)
686     enddo
687     enddo
688     enddo
689     enddo
690     #endif
691    
692     c-- Read meridional wind stress variance.
693 heimbach 1.2 #if (defined (ALLOW_VSTRESS_COST_CONTRIBUTION) || defined (ALLOW_VSTRESS_CONTROL))
694 heimbach 1.1 nnz = 1
695     ce irec = 2
696     ce( due to Patrick's processing:
697     irec = 1
698     ce)
699    
700     if ( tauv_errfile .NE. ' ' ) then
701     call mdsreadfield( tauv_errfile, cost_iprec, cost_yftype, nnz,
702 heimbach 1.2 & wtauv, irec, mythid )
703 heimbach 1.1 endif
704    
705     do bj = jtlo,jthi
706     do bi = itlo,ithi
707     do j = jmin,jmax
708     do i = imin,imax
709     c-- Test for missing values.
710     if (wtauv(i,j,bi,bj) .lt. -9900.) then
711     wtauv(i,j,bi,bj) = 0. _d 0
712     endif
713     wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*0.1
714     wtauv(i,j,bi,bj) = max(wtauv(i,j,bi,bj),wtau0)
715 heimbach 1.2 wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*masks(i,j,k,bi,bj)
716 heimbach 1.1 & *frame(i,j)
717     cph(
718 heimbach 1.2 wtauv2(i,j,bi,bj) = wtau0*maskS(i,j,k,bi,bj)*frame(i,j)
719 heimbach 1.1 cph)
720     enddo
721     enddo
722     enddo
723     enddo
724    
725 heimbach 1.2 #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION) || defined (ALLOW_VWIND_CONTROL))
726 heimbach 1.1
727     nnz = 1
728     ce irec = 2
729     ce( due to Patrick's processing:
730     irec = 1
731     ce)
732    
733     if ( vwind_errfile .NE. ' ' ) then
734     call mdsreadfield( vwind_errfile, cost_iprec, cost_yftype, nnz,
735 heimbach 1.2 & wvwind, irec, mythid )
736 heimbach 1.1 endif
737    
738     do bj = jtlo,jthi
739     do bi = itlo,ithi
740     do j = jmin,jmax
741     do i = imin,imax
742     c-- Test for missing values.
743     if (wvwind(i,j,bi,bj) .lt. -9900.) then
744     wvwind(i,j,bi,bj) = 0. _d 0
745     endif
746     wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)
747     wvwind(i,j,bi,bj) = max(wvwind(i,j,bi,bj),wwind0)
748 heimbach 1.2 wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)*maskc(i,j,k,bi,bj)
749 heimbach 1.1 & *frame(i,j)
750 heimbach 1.2 enddo
751 heimbach 1.1 enddo
752     enddo
753     enddo
754     #endif
755    
756 heimbach 1.2 #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION) || defined (ALLOW_HFLUX_CONTROL))
757 heimbach 1.1 c-- Read heat flux flux variance.
758     nnz = 1
759     c-- First record in data file: mean field.
760     c-- Second record in data file: rms field.
761     ce irec = 2
762     ce( due to Patrick's processing:
763     irec = 1
764     ce)
765     if ( hflux_errfile .NE. ' ' ) then
766 heimbach 1.2 call mdsreadfield( hflux_errfile, cost_iprec, cost_yftype, nnz,
767     & whflux, irec, mythid )
768 heimbach 1.1 endif
769    
770     do bj = jtlo,jthi
771     do bi = itlo,ithi
772     do j = jmin,jmax
773     do i = imin,imax
774     c-- Test for missing values.
775     if (whflux(i,j,bi,bj) .lt. -9900.) then
776     whflux(i,j,bi,bj) = 0. _d 0
777     endif
778     c-- Data are in units of W/m**2.
779     whflux(i,j,bi,bj) = whflux(i,j,bi,bj)/3.
780     whflux(i,j,bi,bj) = max(whflux(i,j,bi,bj),whflux0)
781     & *frame(i,j)
782     whfluxm(i,j,bi,bj) = max(whfluxm(i,j,bi,bj),whflux0m)
783     & *frame(i,j)
784     cph(
785 heimbach 1.2 whflux2(i,j,bi,bj) = whflux0*frame(i,j)
786 heimbach 1.1 cph)
787     enddo
788     enddo
789     enddo
790     enddo
791 heimbach 1.2 #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION) || defined (ALLOW_ATEMP_CONTROL))
792 heimbach 1.1 c-- Read atmos. temp. variance.
793     nnz = 1
794     c-- First record in data file: mean field.
795     c-- Second record in data file: rms field.
796     ce irec = 2
797     ce( due to Patrick's processing:
798     irec = 1
799     ce)
800     if ( atemp_errfile .NE. ' ' ) then
801     call mdsreadfield( atemp_errfile, cost_iprec, cost_yftype, nnz,
802 heimbach 1.2 & watemp, irec, mythid )
803 heimbach 1.1 endif
804    
805     do bj = jtlo,jthi
806     do bi = itlo,ithi
807     do j = jmin,jmax
808     do i = imin,imax
809     c-- Test for missing values.
810     if (watemp(i,j,bi,bj) .lt. -9900.) then
811     watemp(i,j,bi,bj) = 0. _d 0
812     endif
813 heimbach 1.2 c-- Data are in units of W/m**2.
814 heimbach 1.1 watemp(i,j,bi,bj) = watemp(i,j,bi,bj)
815     watemp(i,j,bi,bj) = max(watemp(i,j,bi,bj),watemp0)
816     & *frame(i,j)
817     enddo
818     enddo
819     enddo
820     enddo
821 heimbach 1.2
822    
823 heimbach 1.1 #endif
824    
825 heimbach 1.2 #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SFLUX_CONTROL))
826 heimbach 1.1 c-- Read salt flux variance. Second read: data in units of m/s.
827     nnz = 1
828     c-- First record in data file: mean field.
829     c-- Second record in data file: rms field.
830     ce irec = 2
831     ce( due to Patrick's processing:
832     irec = 1
833     ce)
834     if ( sflux_errfile .NE. ' ' ) then
835 heimbach 1.2 call mdsreadfield( sflux_errfile, cost_iprec, cost_yftype, nnz,
836     & wsflux, irec, mythid )
837 heimbach 1.1 endif
838    
839     do bj = jtlo,jthi
840     do bi = itlo,ithi
841     do j = jmin,jmax
842     do i = imin,imax
843     c-- Test for missing values.
844     if (wsflux(i,j,bi,bj) .lt. -9900.) then
845     wsflux(i,j,bi,bj) = 0. _d 0
846     endif
847     c-- Data are in units of m/s.
848     wsflux(i,j,bi,bj) = wsflux(i,j,bi,bj) / 3.
849     wsflux(i,j,bi,bj) = max(wsflux(i,j,bi,bj),wsflux0)
850     & *frame(i,j)
851     wsfluxm(i,j,bi,bj) = max(wsfluxm(i,j,bi,bj),wsflux0m)
852     & *frame(i,j)
853     cph(
854 heimbach 1.2 wsflux2(i,j,bi,bj) = wsflux0*frame(i,j)
855 heimbach 1.1 cph)
856     enddo
857     enddo
858     enddo
859     enddo
860 heimbach 1.2 #elif (defined (ALLOW_AQH_COST_CONTRIBUTION) || defined (ALLOW_AQH_CONTROL))
861 heimbach 1.1 c-- Secific humid. variance. Second read: data in units of m/s.
862     nnz = 1
863     c-- First record in data file: mean field.
864     c-- Second record in data file: rms field.
865     ce irec = 2
866     ce( due to Patrick's processing:
867     irec = 1
868     ce)
869     if ( aqh_errfile .NE. ' ' ) then
870 heimbach 1.2 call mdsreadfield( aqh_errfile, cost_iprec, cost_yftype, nnz,
871     & waqh, irec, mythid )
872 heimbach 1.1 endif
873    
874     do bj = jtlo,jthi
875     do bi = itlo,ithi
876     do j = jmin,jmax
877     do i = imin,imax
878     c-- Test for missing values.
879     if (waqh(i,j,bi,bj) .lt. -9900.) then
880     waqh(i,j,bi,bj) = 0. _d 0
881     endif
882 heimbach 1.2 c-- Data are in units of m/s.
883 heimbach 1.1 waqh(i,j,bi,bj) = waqh(i,j,bi,bj)
884     waqh(i,j,bi,bj) = max(waqh(i,j,bi,bj),waqh0)
885     & *frame(i,j)
886     enddo
887     enddo
888     enddo
889     enddo
890     #endif
891    
892 heimbach 1.2 #if (defined (ALLOW_PRECIP_COST_CONTRIBUTION) || defined (ALLOW_PRECIP_CONTROL))
893     c-- Atmos. precipitation variance. Second read: data in units of m/s.
894     nnz = 1
895     c-- First record in data file: mean field.
896     c-- Second record in data file: rms field.
897     ce irec = 2
898     ce( due to Patrick's processing:
899     irec = 1
900     ce)
901     if ( precip_errfile .NE. ' ' ) then
902     call mdsreadfield( precip_errfile, cost_iprec, cost_yftype, nnz,
903     & wprecip, irec, mythid )
904     endif
905    
906     do bj = jtlo,jthi
907     do bi = itlo,ithi
908     do j = jmin,jmax
909     do i = imin,imax
910     c-- Test for missing values.
911     if (wprecip(i,j,bi,bj) .lt. -9900.) then
912     wprecip(i,j,bi,bj) = 0. _d 0
913     endif
914     c-- Data are in units of m/s.
915     wprecip(i,j,bi,bj) = wprecip(i,j,bi,bj)
916     wprecip(i,j,bi,bj) = max(wprecip(i,j,bi,bj),wprecip0)
917     & *frame(i,j)
918     enddo
919     enddo
920     enddo
921     enddo
922     #endif
923    
924     #if (defined (ALLOW_SWFLUX_COST_CONTRIBUTION) || defined (ALLOW_SWFLUX_CONTROL))
925     c-- Atmos. swfluxitation variance. Second read: data in units of m/s.
926     nnz = 1
927     c-- First record in data file: mean field.
928     c-- Second record in data file: rms field.
929     ce irec = 2
930     ce( due to Patrick's processing:
931     irec = 1
932     ce)
933     if ( swflux_errfile .NE. ' ' ) then
934     call mdsreadfield( swflux_errfile, cost_iprec, cost_yftype, nnz,
935     & wswflux, irec, mythid )
936     endif
937    
938     do bj = jtlo,jthi
939     do bi = itlo,ithi
940     do j = jmin,jmax
941     do i = imin,imax
942     c-- Test for missing values.
943     if (wswflux(i,j,bi,bj) .lt. -9900.) then
944     wswflux(i,j,bi,bj) = 0. _d 0
945     endif
946     c-- Data are in units of m/s.
947     wswflux(i,j,bi,bj) = wswflux(i,j,bi,bj)
948     wswflux(i,j,bi,bj) = max(wswflux(i,j,bi,bj),wswflux0)
949     & *frame(i,j)
950     enddo
951     enddo
952     enddo
953     enddo
954     #endif
955    
956     #if (defined (ALLOW_SWDOWN_COST_CONTRIBUTION) || defined (ALLOW_SWDOWN_CONTROL))
957     c-- Atmos. swdownitation variance. Second read: data in units of m/s.
958     nnz = 1
959     c-- First record in data file: mean field.
960     c-- Second record in data file: rms field.
961     ce irec = 2
962     ce( due to Patrick's processing:
963     irec = 1
964     ce)
965     if ( swdown_errfile .NE. ' ' ) then
966     call mdsreadfield( swdown_errfile, cost_iprec, cost_yftype, nnz,
967     & wswdown, irec, mythid )
968     endif
969    
970     do bj = jtlo,jthi
971     do bi = itlo,ithi
972     do j = jmin,jmax
973     do i = imin,imax
974     c-- Test for missing values.
975     if (wswdown(i,j,bi,bj) .lt. -9900.) then
976     wswdown(i,j,bi,bj) = 0. _d 0
977     endif
978     c-- Data are in units of m/s.
979     wswdown(i,j,bi,bj) = wswdown(i,j,bi,bj)
980     wswdown(i,j,bi,bj) = max(wswdown(i,j,bi,bj),wswdown0)
981     & *frame(i,j)
982     enddo
983     enddo
984     enddo
985     enddo
986     #endif
987    
988 heimbach 1.1 c-- Units have to be checked!
989     do bj = jtlo,jthi
990     do bi = itlo,ithi
991     do j = jmin,jmax
992     do i = imin,imax
993     if (wtp(i,j,bi,bj) .ne. 0.) then
994     wtp (i,j,bi,bj) = 1./wtp(i,j,bi,bj)/wtp(i,j,bi,bj)
995     endif
996     if (wers(i,j,bi,bj) .ne. 0.) then
997     wers(i,j,bi,bj) = 1./wers(i,j,bi,bj)/wers(i,j,bi,bj)
998     endif
999     if (wsst(i,j,bi,bj) .ne. 0.) then
1000 heimbach 1.2 wsst(i,j,bi,bj) = 1./wsst(i,j,bi,bj)/wsst(i,j,bi,bj)
1001 heimbach 1.1 endif
1002 heimbach 1.2 if (wsss(i,j,bi,bj) .ne. 0.) then
1003     wsss(i,j,bi,bj) = 1./wsss(i,j,bi,bj)/wsss(i,j,bi,bj)
1004 heimbach 1.1 endif
1005     if (wp(i,j,bi,bj) .ne. 0.) then
1006     wp(i,j,bi,bj) = 1./wp(i,j,bi,bj)/wp(i,j,bi,bj)
1007     endif
1008     if (wtauu(i,j,bi,bj) .ne. 0.) then
1009 heimbach 1.2 wtauu(i,j,bi,bj) = 1./wtauu(i,j,bi,bj)/wtauu(i,j,bi,bj)
1010 heimbach 1.1 else
1011     wtauu(i,j,bi,bj) = 0.0 _d 0
1012     endif
1013     if (wtauum(i,j,bi,bj) .ne. 0.) then
1014     wtauum(i,j,bi,bj) =
1015     & 1./wtauum(i,j,bi,bj)/wtauum(i,j,bi,bj)
1016     else
1017     wtauum(i,j,bi,bj) = 0.0 _d 0
1018     endif
1019     if (wscatx(i,j,bi,bj) .ne. 0.) then
1020     wscatx(i,j,bi,bj) =
1021     & 1./wscatx(i,j,bi,bj)/wscatx(i,j,bi,bj)
1022     else
1023     wscatx(i,j,bi,bj) = 0.0 _d 0
1024     endif
1025     if (wtauv(i,j,bi,bj) .ne. 0.) then
1026 heimbach 1.2 wtauv(i,j,bi,bj) = 1./wtauv(i,j,bi,bj)/wtauv(i,j,bi,bj)
1027 heimbach 1.1 else
1028     wtauv(i,j,bi,bj) = 0.0 _d 0
1029     endif
1030     if (wtauvm(i,j,bi,bj) .ne. 0.) then
1031     wtauvm(i,j,bi,bj) =
1032     & 1./wtauvm(i,j,bi,bj)/wtauvm(i,j,bi,bj)
1033     else
1034     wtauvm(i,j,bi,bj) = 0.0 _d 0
1035     endif
1036     if (wscaty(i,j,bi,bj) .ne. 0.) then
1037     wscaty(i,j,bi,bj) =
1038     & 1./wscaty(i,j,bi,bj)/wscaty(i,j,bi,bj)
1039     else
1040     wscaty(i,j,bi,bj) = 0.0 _d 0
1041     endif
1042     if (whflux(i,j,bi,bj) .ne. 0.) then
1043     whflux(i,j,bi,bj) =
1044     & 1./whflux(i,j,bi,bj)/whflux(i,j,bi,bj)
1045     else
1046     whflux(i,j,bi,bj) = 0.0 _d 0
1047     endif
1048     if (whfluxm(i,j,bi,bj) .ne. 0.) then
1049     whfluxm(i,j,bi,bj) =
1050     & 1./whfluxm(i,j,bi,bj)/whfluxm(i,j,bi,bj)
1051     else
1052     whfluxm(i,j,bi,bj) = 0.0 _d 0
1053     endif
1054     if (wsflux(i,j,bi,bj) .ne. 0.) then
1055     wsflux(i,j,bi,bj) =
1056     & 1./wsflux(i,j,bi,bj)/wsflux(i,j,bi,bj)
1057     else
1058     wsflux(i,j,bi,bj) = 0.0 _d 0
1059     endif
1060     if (wsfluxm(i,j,bi,bj) .ne. 0.) then
1061     wsfluxm(i,j,bi,bj) =
1062     & 1./wsfluxm(i,j,bi,bj)/wsfluxm(i,j,bi,bj)
1063     else
1064     wsfluxm(i,j,bi,bj) = 0.0 _d 0
1065     endif
1066     if (wuwind(i,j,bi,bj) .ne. 0.) then
1067     wuwind(i,j,bi,bj) =
1068     & 1./wuwind(i,j,bi,bj)/wuwind(i,j,bi,bj)
1069     else
1070     wuwind(i,j,bi,bj) = 0.0 _d 0
1071     endif
1072     if (wvwind(i,j,bi,bj) .ne. 0.) then
1073     wvwind(i,j,bi,bj) =
1074     & 1./wvwind(i,j,bi,bj)/wvwind(i,j,bi,bj)
1075     else
1076     wvwind(i,j,bi,bj) = 0.0 _d 0
1077     endif
1078     if (watemp(i,j,bi,bj) .ne. 0.) then
1079     watemp(i,j,bi,bj) =
1080     & 1./watemp(i,j,bi,bj)/watemp(i,j,bi,bj)
1081     else
1082     watemp(i,j,bi,bj) = 0.0 _d 0
1083     endif
1084     if (waqh(i,j,bi,bj) .ne. 0.) then
1085     waqh(i,j,bi,bj) =
1086     & 1./waqh(i,j,bi,bj)/waqh(i,j,bi,bj)
1087     else
1088     waqh(i,j,bi,bj) = 0.0 _d 0
1089     endif
1090 heimbach 1.2 if (wprecip(i,j,bi,bj) .ne. 0.) then
1091     wprecip(i,j,bi,bj) =
1092     & 1./wprecip(i,j,bi,bj)/wprecip(i,j,bi,bj)
1093     else
1094     wprecip(i,j,bi,bj) = 0.0 _d 0
1095     endif
1096     if (wswflux(i,j,bi,bj) .ne. 0.) then
1097     wswflux(i,j,bi,bj) =
1098     & 1./wswflux(i,j,bi,bj)/wswflux(i,j,bi,bj)
1099     else
1100     wswflux(i,j,bi,bj) = 0.0 _d 0
1101     endif
1102     if (wswdown(i,j,bi,bj) .ne. 0.) then
1103     wswdown(i,j,bi,bj) =
1104     & 1./wswdown(i,j,bi,bj)/wswdown(i,j,bi,bj)
1105     else
1106     wswdown(i,j,bi,bj) = 0.0 _d 0
1107     endif
1108 heimbach 1.1
1109     if (wsfluxmm(bi,bj).ne.0.)
1110     & wsfluxmm(bi,bj) = 1./wsfluxmm(bi,bj)*wsfluxmm(bi,bj)
1111     if (whfluxmm(bi,bj).ne.0.)
1112     & whfluxmm(bi,bj) = 1./whfluxmm(bi,bj)*whfluxmm(bi,bj)
1113 heimbach 1.2
1114 heimbach 1.1 cph(
1115     if (whflux2(i,j,bi,bj) .ne. 0.) then
1116     whflux2(i,j,bi,bj) =
1117     & 1./whflux2(i,j,bi,bj)/whflux2(i,j,bi,bj)
1118     else
1119     whflux2(i,j,bi,bj) = 0.0 _d 0
1120     endif
1121     if (wsflux2(i,j,bi,bj) .ne. 0.) then
1122     wsflux2(i,j,bi,bj) =
1123     & 1./wsflux2(i,j,bi,bj)/wsflux2(i,j,bi,bj)
1124     else
1125     wsflux2(i,j,bi,bj) = 0.0 _d 0
1126     endif
1127     if (wtauu2(i,j,bi,bj) .ne. 0.) then
1128 heimbach 1.2 wtauu2(i,j,bi,bj) =
1129 heimbach 1.1 & 1./wtauu2(i,j,bi,bj)/wtauu2(i,j,bi,bj)
1130     else
1131     wtauu2(i,j,bi,bj) = 0.0 _d 0
1132     endif
1133     if (wtauv2(i,j,bi,bj) .ne. 0.) then
1134     wtauv2(i,j,bi,bj) =
1135     & 1./wtauv2(i,j,bi,bj)/wtauv2(i,j,bi,bj)
1136     else
1137     wtauv2(i,j,bi,bj) = 0.0 _d 0
1138     endif
1139     cph)
1140     enddo
1141     enddo
1142    
1143     do k = 1,nr
1144     if (wtheta(k,bi,bj) .ne. 0.) then
1145 heimbach 1.2 wtheta(k,bi,bj) = ratio/wtheta(k,bi,bj)/wtheta(k,bi,bj)
1146 heimbach 1.1 else
1147     wtheta(k,bi,bj) = 0.0 _d 0
1148     endif
1149     if (wsalt(k,bi,bj) .ne. 0.) then
1150 heimbach 1.2 wsalt(k,bi,bj) = ratio/wsalt(k,bi,bj)/wsalt(k,bi,bj)
1151 heimbach 1.1 else
1152     wsalt(k,bi,bj) = 0.0 _d 0
1153     endif
1154     enddo
1155    
1156     #ifdef ALLOW_OBCS_COST_CONTRIBUTION
1157     do iobcs = 1,nobcs
1158     do k = 1,nr
1159     #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
1160     if (wobcsn(k,iobcs) .ne. 0.) then
1161     wobcsn(k,iobcs) =
1162     & ratio/wobcsn(k,iobcs)/wobcsn(k,iobcs)
1163     else
1164     wobcsn(k,iobcs) = 0.0 _d 0
1165     endif
1166     #endif
1167     #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
1168     if (wobcss(k,iobcs) .ne. 0.) then
1169     wobcss(k,iobcs) =
1170     & ratio/wobcss(k,iobcs)/wobcss(k,iobcs)
1171     else
1172     wobcss(k,iobcs) = 0.0 _d 0
1173     endif
1174     #endif
1175     #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
1176     if (wobcsw(k,iobcs) .ne. 0.) then
1177     wobcsw(k,iobcs) =
1178     & ratio/wobcsw(k,iobcs)/wobcsw(k,iobcs)
1179     else
1180     wobcsw(k,iobcs) = 0.0 _d 0
1181     endif
1182     #endif
1183     #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
1184     if (wobcse(k,iobcs) .ne. 0.) then
1185     wobcse(k,iobcs) =
1186     & ratio/wobcse(k,iobcs)/wobcse(k,iobcs)
1187     else
1188     wobcse(k,iobcs) = 0.0 _d 0
1189     endif
1190     #endif
1191     enddo
1192     enddo
1193     #endif
1194    
1195     enddo
1196     enddo
1197    
1198 heimbach 1.2 do bj = jtlo,jthi
1199     do bi = itlo,ithi
1200     do k = 1,nr
1201     if (wdiffkr(k,bi,bj) .ne. 0.) then
1202     wdiffkr(k,bi,bj) = 1./wdiffkr(k,bi,bj)/wdiffkr(k,bi,bj)
1203     else
1204     wdiffkr(k,bi,bj) = 0.0 _d 0
1205     endif
1206     if (wkapgm(k,bi,bj) .ne. 0.) then
1207     wkapgm(k,bi,bj) = 1./wkapgm(k,bi,bj)/wkapgm(k,bi,bj)
1208     else
1209     wkapgm(k,bi,bj) = 0.0 _d 0
1210     endif
1211     if (wedtaux(k,bi,bj) .ne. 0.) then
1212     wedtaux(k,bi,bj) = 1./wedtaux(k,bi,bj)/wedtaux(k,bi,bj)
1213     else
1214     wedtaux(k,bi,bj) = 0.0 _d 0
1215     endif
1216     if (wedtauy(k,bi,bj) .ne. 0.) then
1217     wedtauy(k,bi,bj) = 1./wedtauy(k,bi,bj)/wedtauy(k,bi,bj)
1218     else
1219     wedtauy(k,bi,bj) = 0.0 _d 0
1220     endif
1221     do j = jmin,jmax
1222     do i = imin,imax
1223     if (wdiffkr2(i,j,k,bi,bj) .ne. 0.) then
1224     wdiffkr2(i,j,k,bi,bj) = frame(i,j)/
1225     & wdiffkr2(i,j,k,bi,bj)/wdiffkr2(i,j,k,bi,bj)
1226     else
1227     wdiffkr2(i,j,k,bi,bj) = 0.0 _d 0
1228     endif
1229     wdiffkrFld(i,j,k,bi,bj) = wdiffkr2(i,j,k,bi,bj)
1230     c
1231     if (wkapgm2(i,j,k,bi,bj) .ne. 0.) then
1232     wkapgm2(i,j,k,bi,bj) = frame(i,j)/
1233     & wkapgm2(i,j,k,bi,bj)/wkapgm2(i,j,k,bi,bj)
1234     else
1235     wkapgm2(i,j,k,bi,bj) = 0.0 _d 0
1236     endif
1237     wkapgmFld(i,j,k,bi,bj) = wkapgm2(i,j,k,bi,bj)
1238     c
1239     if (wedtaux2(i,j,k,bi,bj) .ne. 0.) then
1240     wedtaux2(i,j,k,bi,bj) = frame(i,j)/
1241     & wedtaux2(i,j,k,bi,bj)/wedtaux2(i,j,k,bi,bj)
1242     else
1243     wedtaux2(i,j,k,bi,bj) = 0.0 _d 0
1244     endif
1245     wedtauxFld(i,j,k,bi,bj) = wedtaux2(i,j,k,bi,bj)
1246     c
1247     if (wedtauy2(i,j,k,bi,bj) .ne. 0.) then
1248     wedtauy2(i,j,k,bi,bj) = frame(i,j)/
1249     & wedtauy2(i,j,k,bi,bj)/wedtauy2(i,j,k,bi,bj)
1250     else
1251     wedtauy2(i,j,k,bi,bj) = 0.0 _d 0
1252     endif
1253     wedtauyFld(i,j,k,bi,bj) = wedtauy2(i,j,k,bi,bj)
1254     enddo
1255     enddo
1256     enddo
1257     enddo
1258     enddo
1259     c
1260     c write masks and weights to files to be read by a master process
1261     c
1262     call active_write_xyz_loc( 'maskCtrlC', maskC,
1263     & 1, 0, mythid, dummy)
1264     call active_write_xyz_loc( 'maskCtrlW', maskW,
1265     & 1, 0, mythid, dummy)
1266     call active_write_xyz_loc( 'maskCtrlS', maskS,
1267     & 1, 0, mythid, dummy)
1268    
1269     #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION) || defined (ALLOW_HFLUX_CONTROL))
1270 heimbach 1.1 call active_write_xy_loc( 'whflux', whflux, 1, 0, mythid, dummy)
1271     call active_write_xy_loc( 'whflux2', whflux2, 1, 0, mythid, dummy)
1272 heimbach 1.2 #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION) || defined (ALLOW_ATEMP_CONTROL))
1273 heimbach 1.1 call active_write_xy_loc( 'watemp', watemp, 1, 0, mythid, dummy)
1274     #endif
1275    
1276 heimbach 1.2 #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SFLUX_CONTROL))
1277 heimbach 1.1 call active_write_xy_loc( 'wsflux', wsflux, 1, 0, mythid, dummy)
1278     call active_write_xy_loc( 'wsflux2', wsflux2, 1, 0, mythid, dummy)
1279 heimbach 1.2 #elif (defined (ALLOW_AQH_COST_CONTRIBUTION) || defined (ALLOW_AQH_CONTROL))
1280 heimbach 1.1 call active_write_xy_loc( 'waqh', waqh, 1, 0, mythid, dummy)
1281     #endif
1282    
1283 heimbach 1.2 #if (defined (ALLOW_PRECIP_COST_CONTRIBUTION) || defined (ALLOW_PRECIP_CONTROL))
1284     call active_write_xy_loc( 'wprecip', wprecip, 1, 0, mythid, dummy)
1285     #endif
1286    
1287     #if (defined (ALLOW_SWFLUX_COST_CONTRIBUTION) || defined (ALLOW_SWFLUX_CONTROL))
1288     call active_write_xy_loc( 'wswflux', wswflux, 1, 0, mythid, dummy)
1289     #endif
1290    
1291     #if (defined (ALLOW_SWDOWN_COST_CONTRIBUTION) || defined (ALLOW_SWDOWN_CONTROL))
1292     call active_write_xy_loc( 'wswdown', wswdown, 1, 0, mythid, dummy)
1293     #endif
1294    
1295     #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION) || defined (ALLOW_USTRESS_CONTROL))
1296 heimbach 1.1 call active_write_xy_loc( 'wtauu', wtauu, 1, 0, mythid, dummy)
1297     call active_write_xy_loc( 'wtauu2', wtauu2, 1, 0, mythid, dummy)
1298 heimbach 1.2 #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION) || defined (ALLOW_UWIND_CONTROL))
1299 heimbach 1.1 call active_write_xy_loc( 'wuwind', wuwind, 1, 0, mythid, dummy)
1300     #endif
1301    
1302 heimbach 1.2 #if (defined (ALLOW_VSTRESS_COST_CONTRIBUTION) || defined (ALLOW_VSTRESS_CONTROL))
1303 heimbach 1.1 call active_write_xy_loc( 'wtauv', wtauv, 1, 0, mythid, dummy)
1304     call active_write_xy_loc( 'wtauv2', wtauv2, 1, 0, mythid, dummy)
1305 heimbach 1.2 #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION) || defined (ALLOW_VWIND_CONTROL))
1306 heimbach 1.1 call active_write_xy_loc( 'wvwind', wvwind, 1, 0, mythid, dummy)
1307     #endif
1308    
1309 heimbach 1.2 #if (defined (ALLOW_SST_COST_CONTRIBUTION) || defined (ALLOW_SST_CONTROL))
1310     call active_write_xy_loc( 'wsst', wsst, 1, 0, mythid, dummy)
1311     #endif
1312    
1313     #if (defined (ALLOW_SSS_COST_CONTRIBUTION) || defined (ALLOW_SSS_CONTROL))
1314     call active_write_xy_loc( 'wsss', wsss, 1, 0, mythid, dummy)
1315     #endif
1316    
1317 heimbach 1.1 #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
1318     #endif
1319    
1320     end
  ViewVC Help
Powered by ViewVC 1.1.22