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

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

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

revision 1.1 by heimbach, Wed May 3 23:59:39 2006 UTC revision 1.2 by heimbach, Wed May 31 21:30:36 2006 UTC
# Line 86  c     == end of interface == Line 86  c     == end of interface ==
86        imax = snx+olx        imax = snx+olx
87    
88  c--   Initialize background weights  c--   Initialize background weights
89        wtau0    = 0.        whflux0m = whflux0
90        whflux0  = 0.        wsflux0m = wsflux0
91        wsflux0  = 0.        wtau0m   = wtau0
       whflux0m = 0  
       wsflux0m = 0.  
       watemp0  = 0.  
       waqh0    = 0.  
       wwind0   = 0.  
92    
93  c--   Initialize variance (weight) fields.  c--   Initialize variance (weight) fields.
94        do k = 1,nr        do k = 1,nr
# Line 116  c--   Initialize variance (weight) field Line 111  c--   Initialize variance (weight) field
111                wtauvm  (i,j,bi,bj) = 0. _d 0                wtauvm  (i,j,bi,bj) = 0. _d 0
112                watemp  (i,j,bi,bj) = 0. _d 0                watemp  (i,j,bi,bj) = 0. _d 0
113                waqh    (i,j,bi,bj) = 0. _d 0                waqh    (i,j,bi,bj) = 0. _d 0
114                  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                wuwind  (i,j,bi,bj) = 0. _d 0                wuwind  (i,j,bi,bj) = 0. _d 0
118                wvwind  (i,j,bi,bj) = 0. _d 0                wvwind  (i,j,bi,bj) = 0. _d 0
119                wsst    (i,j,bi,bj) = 0. _d 0                wsst    (i,j,bi,bj) = 0. _d 0
120                wsss    (i,j,bi,bj) = 0. _d 0                wsss    (i,j,bi,bj) = 0. _d 0
121                wtp     (i,j,bi,bj) = 0. _d 0                wtp     (i,j,bi,bj) = 0. _d 0
122                wers    (i,j,bi,bj) = 0. _d 0                wers    (i,j,bi,bj) = 0. _d 0
               wgfo    (i,j,bi,bj) = 0. _d 0  
123                wp      (i,j,bi,bj) = 0. _d 0                wp      (i,j,bi,bj) = 0. _d 0
124                wudrift (i,j,bi,bj) = 0. _d 0                wudrift (i,j,bi,bj) = 0. _d 0
125                wvdrift (i,j,bi,bj) = 0. _d 0                wvdrift (i,j,bi,bj) = 0. _d 0
126  cph(  cph(
127                whflux2 (i,j,bi,bj) = 0. _d 0                whflux2 (i,j,bi,bj) = 0. _d 0
128                wsflux2 (i,j,bi,bj) = 0. _d 0                wsflux2 (i,j,bi,bj) = 0. _d 0
129                wtauu2  (i,j,bi,bj) = 0. _d 0                wtauu2  (i,j,bi,bj) = 0. _d 0
130                wtauv2  (i,j,bi,bj) = 0. _d 0                wtauv2  (i,j,bi,bj) = 0. _d 0
131  cph)  cph)
132              enddo              enddo
133            enddo            enddo
# Line 143  cph) Line 140  cph)
140              wsalt  (k,bi,bj) = 0. _d 0              wsalt  (k,bi,bj) = 0. _d 0
141              wctdt  (k,bi,bj) = 0. _d 0              wctdt  (k,bi,bj) = 0. _d 0
142              wctds  (k,bi,bj) = 0. _d 0              wctds  (k,bi,bj) = 0. _d 0
143                wdiffkr(k,bi,bj) = wdiffkr0
144                wkapgm (k,bi,bj) = wkapgm0
145                wedtaux(k,bi,bj) = wedtau0
146                wedtauy(k,bi,bj) = wedtau0
147              do j = jmin,jmax              do j = jmin,jmax
148                do i = imin,imax                do i = imin,imax
149                  wtheta2 (i,j,k,bi,bj) = 0. _d 0                  wtheta2 (i,j,k,bi,bj) = 0. _d 0
150                  wsalt2  (i,j,k,bi,bj) = 0. _d 0                  wsalt2  (i,j,k,bi,bj) = 0. _d 0
151                    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                  wthetaLev (i,j,k,bi,bj) = 0. _d 0                  wthetaLev (i,j,k,bi,bj) = 0. _d 0
156                  wsaltLev  (i,j,k,bi,bj) = 0. _d 0                  wsaltLev  (i,j,k,bi,bj) = 0. _d 0
157                    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                enddo                enddo
162              enddo              enddo
163            enddo            enddo
# Line 218  c--   Read error information and set up Line 227  c--   Read error information and set up
227        _BEGIN_MASTER(myThid)        _BEGIN_MASTER(myThid)
228          ilo = ifnblnk(data_errfile)          ilo = ifnblnk(data_errfile)
229          ihi = ilnblnk(data_errfile)          ihi = ilnblnk(data_errfile)
230          CALL OPEN_COPY_DATA_FILE(          CALL OPEN_COPY_DATA_FILE(
231       I                          data_errfile(ilo:ihi),       I                          data_errfile(ilo:ihi),
232       I                          'ECCO_COST_WEIGHTS',       I                          'ECCO_COST_WEIGHTS',
233       O                          gwunit,       O                          gwunit,
234       I                          myThid )       I                          myThid )
235    
236          read(gwunit,*)          read(gwunit,*) ratio
 #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION) || defined (ALLOW_HFLUX_CONTROL))  
      &         whflux0  
 #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION) || defined (ALLOW_ATEMP_CONTROL))  
      &         watemp0  
 #endif  
 #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SFLUX_CONTROL))  
      &       , wsflux0  
 #elif (defined (ALLOW_AQH_COST_CONTRIBUTION) || defined (ALLOW_AQH_CONTROL))  
      &       , waqh0  
 #endif  
 #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION) || defined (ALLOW_USTRESS_CONTROL))  
      &       , wtau0  
 #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION) || defined (ALLOW_UWIND_CONTROL))  
      &       , wwind0  
 #endif  
      &       , ratio  
237  #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || defined (ALLOW_OBCS_CONTROL))  #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || defined (ALLOW_OBCS_CONTROL))
238       &       , wbaro       &       , wbaro
239  #endif  #endif
   
240          do k = 1,nr          do k = 1,nr
241            read(gwunit,*) wti(k), wsi(k)            read(gwunit,*) wti(k), wsi(k)
242  #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || \  #if (defined (ALLOW_OBCS_COST_CONTRIBUTION) || defined (ALLOW_OBCS_CONTROL))
      defined (ALLOW_OBCS_CONTROL))  
243       &               , wvi(k)       &               , wvi(k)
244  #endif  #endif
245          end do          end do
246          close(gwunit)          close(gwunit)
247    
         whflux0m = whflux0  
         wsflux0m = wsflux0  
         wtau0m   = wtau0  
248        _END_MASTER(myThid)        _END_MASTER(myThid)
249    
250        _BARRIER        _BARRIER
# Line 285  c--       is no need to add the zonal an Line 273  c--       is no need to add the zonal an
273  #ifdef ALLOW_OBCSN_COST_CONTRIBUTION  #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
274              wobcsn(k,1) = wti(k)              wobcsn(k,1) = wti(k)
275              wobcsn(k,2) = wsi(k)              wobcsn(k,2) = wsi(k)
276              wobcsn(k,3) = wti(k)*0.02              wobcsn(k,3) = wvi(k)
277              wobcsn(k,4) = wti(k)*0.02              wobcsn(k,4) = wvi(k)
278  #endif  #endif
279  #ifdef ALLOW_OBCSS_COST_CONTRIBUTION  #ifdef ALLOW_OBCSS_COST_CONTRIBUTION
280              wobcss(k,1) = wti(k)              wobcss(k,1) = wti(k)
281              wobcss(k,2) = wsi(k)              wobcss(k,2) = wsi(k)
282              wobcss(k,3) = wti(k)*0.02              wobcss(k,3) = wvi(k)
283              wobcss(k,4) = wti(k)*0.02              wobcss(k,4) = wvi(k)
284  #endif  #endif
285  #ifdef ALLOW_OBCSW_COST_CONTRIBUTION  #ifdef ALLOW_OBCSW_COST_CONTRIBUTION
286              wobcsw(k,1) = wti(k)              wobcsw(k,1) = wti(k)
287              wobcsw(k,2) = wsi(k)              wobcsw(k,2) = wsi(k)
288              wobcsw(k,3) = wti(k)*0.02              wobcsw(k,3) = wvi(k)
289              wobcsw(k,4) = wti(k)*0.02              wobcsw(k,4) = wvi(k)
290  #endif  #endif
291  #ifdef ALLOW_OBCSE_COST_CONTRIBUTION  #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
292              wobcse(k,1) = wti(k)              wobcse(k,1) = wti(k)
293              wobcse(k,2) = wsi(k)              wobcse(k,2) = wsi(k)
294              wobcse(k,3) = wti(k)*0.02              wobcse(k,3) = wvi(k)
295              wobcse(k,4) = wti(k)*0.02              wobcse(k,4) = wvi(k)
296  #endif  #endif
297            enddo            enddo
298          enddo          enddo
# Line 323  c--           Test for missing values. Line 311  c--           Test for missing values.
311                if ( wsaltLev(i,j,k,bi,bj).eq.0 ) then                if ( wsaltLev(i,j,k,bi,bj).eq.0 ) then
312                   wsaltLev(i,j,k,bi,bj) = 0. _d 0                   wsaltLev(i,j,k,bi,bj) = 0. _d 0
313                else                else
314                   wsaltLev(i,j,k,bi,bj)=frame(i,j)/                   wsaltLev(i,j,k,bi,bj)=frame(i,j)/(
315       $              ( wsaltLev(i,j,k,bi,bj)*wsaltLev(i,j,k,bi,bj) )       $                wsaltLev(i,j,k,bi,bj)*wsaltLev(i,j,k,bi,bj) )
316                endif                endif
317               enddo               enddo
318              enddo              enddo
# Line 343  c--           Test for missing values. Line 331  c--           Test for missing values.
331              enddo              enddo
332             enddo             enddo
333            enddo            enddo
334           enddo                       enddo
335        endif        endif
336        call active_write_xyz( 'wsaltLev', wsaltLev,        call active_write_xyz( 'wsaltLev', wsaltLev,
337       &     1, 0, mythid, dummy)       &     1, 0, mythid, dummy)
338  #endif  #endif
339    
# Line 362  c--           Test for missing values. Line 350  c--           Test for missing values.
350                if ( wthetaLev(i,j,k,bi,bj).eq.0 ) then                if ( wthetaLev(i,j,k,bi,bj).eq.0 ) then
351                   wthetaLev(i,j,k,bi,bj) = 0. _d 0                   wthetaLev(i,j,k,bi,bj) = 0. _d 0
352                else                else
353                   wthetaLev(i,j,k,bi,bj)=frame(i,j)/                   wthetaLev(i,j,k,bi,bj)=frame(i,j)/(
354       $              ( wthetaLev(i,j,k,bi,bj)*wthetaLev(i,j,k,bi,bj) )       $                wthetaLev(i,j,k,bi,bj)*wthetaLev(i,j,k,bi,bj) )
355                endif                endif
356               enddo               enddo
357              enddo              enddo
# Line 382  c--           Test for missing values. Line 370  c--           Test for missing values.
370              enddo              enddo
371             enddo             enddo
372            enddo            enddo
373           enddo                       enddo
374        endif        endif
375        call active_write_xyz( 'wthetaLev', wthetaLev,        call active_write_xyz( 'wthetaLev', wthetaLev,
376       &     1, 0, mythid, dummy)       &     1, 0, mythid, dummy)
377  #endif  #endif
378    
# Line 406  c--           Test for missing values. Line 394  c--           Test for missing values.
394       $            wsalt2(i,j,k,bi,bj).eq.0) then       $            wsalt2(i,j,k,bi,bj).eq.0) then
395                   wsalt2(i,j,k,bi,bj) = 0. _d 0                   wsalt2(i,j,k,bi,bj) = 0. _d 0
396                else                else
397  cph  new weights by G. Forget dont need MAX                   wsalt2(i,j,k,bi,bj)=frame(i,j)/
398                   wsalt2(i,j,k,bi,bj)=frame(i,j)/(       $                  ( wsalt2(i,j,k,bi,bj)*wsalt2(i,j,k,bi,bj) )
      $                wsalt2(i,j,k,bi,bj)*wsalt2(i,j,k,bi,bj) )  
399                endif                endif
400               enddo               enddo
401              enddo              enddo
# Line 449  c--           Test for missing values. Line 436  c--           Test for missing values.
436       $            wtheta2(i,j,k,bi,bj).eq.0) then       $            wtheta2(i,j,k,bi,bj).eq.0) then
437                   wtheta2(i,j,k,bi,bj) = 0. _d 0                   wtheta2(i,j,k,bi,bj) = 0. _d 0
438                else                else
439  cph  new weights by G. Forget dont need MAX                   wtheta2(i,j,k,bi,bj)= frame(i,j)/
440                   wtheta2(i,j,k,bi,bj)=frame(i,j)/(       $                  ( wtheta2(i,j,k,bi,bj)*wtheta2(i,j,k,bi,bj) )
      $                wtheta2(i,j,k,bi,bj)*wtheta2(i,j,k,bi,bj) )  
441                endif                endif
442               enddo               enddo
443              enddo              enddo
# Line 510  c--   Set all tile edges to zero. Line 496  c--   Set all tile edges to zero.
496              do i = imin,imax              do i = imin,imax
497                wp(i,j,bi,bj) = wp(i,j,bi,bj)*frame(i,j)                wp(i,j,bi,bj) = wp(i,j,bi,bj)*frame(i,j)
498  cph-indonesian(  cph-indonesian(
499                if ( xC(i,j,bi,bj) .GT. 115. .AND.                if ( xC(i,j,bi,bj) .GT. 120. .AND.
500       &             xC(i,j,bi,bj) .LT. 130. .AND.       &             xC(i,j,bi,bj) .LT. 130. .AND.
501       &             yC(i,j,bi,bj) .GT. -10. .AND.       &             yC(i,j,bi,bj) .GT. -10. .AND.
502       &             yC(i,j,bi,bj) .LT.  10. ) then       &             yC(i,j,bi,bj) .LT.  10. ) then
503  cph                 wp(i,j,bi,bj) = wp(i,j,bi,bj)*10000.                   wp(i,j,bi,bj) = wp(i,j,bi,bj)*100.
                  wp(i,j,bi,bj) = 0.  
504                endif                endif
505  cph-indonesian)  cph-indonesian)
 cph-medit(  
               if ( ( xC(i,j,bi,bj) .GT. 355. .AND.  
      &               xC(i,j,bi,bj) .LT. 360. .AND.  
      &               yC(i,j,bi,bj) .GT.  30. .AND.  
      &               yC(i,j,bi,bj) .LT.  48. )  
      &             .OR.  
      &             ( xC(i,j,bi,bj) .GT.   0. .AND.  
      &               xC(i,j,bi,bj) .LT.  39. .AND.  
      &               yC(i,j,bi,bj) .GT.  30. .AND.  
      &               yC(i,j,bi,bj) .LT.  48. ) ) then  
                  wp(i,j,bi,bj) = wp(i,j,bi,bj)*10.  
               endif  
 cph-medit)  
506              enddo              enddo
507            enddo            enddo
508          enddo          enddo
# Line 564  c--           T/P error + 5cm Line 536  c--           T/P error + 5cm
536                if (_hFacC(i,j,k,bi,bj) .eq. 0.) then                if (_hFacC(i,j,k,bi,bj) .eq. 0.) then
537                  wtp (i,j,bi,bj) = 0. _d 0                  wtp (i,j,bi,bj) = 0. _d 0
538                  wers(i,j,bi,bj) = 0. _d 0                  wers(i,j,bi,bj) = 0. _d 0
                 wgfo(i,j,bi,bj) = 0. _d 0  
539                else                else
540                  wtp (i,j,bi,bj) = ( wtp(i,j,bi,bj) * 0.01 * 0.5 )                  wtp (i,j,bi,bj) = ( wtp(i,j,bi,bj) * 0.01 * 0.5 )
541       &                            *frame(i,j)       &                            *frame(i,j)
542                  wers(i,j,bi,bj) = ( wtp(i,j,bi,bj) + 0.05 )                  wers(i,j,bi,bj) = ( wtp(i,j,bi,bj) + 0.05 )
543       &                            *frame(i,j)       &                            *frame(i,j)
                 wgfo(i,j,bi,bj) = ( wtp(i,j,bi,bj) + 0.05 )  
      &                            *frame(i,j)  
544                endif                endif
545  cph-indonesian(  cph-indonesian(
546                if ( xC(i,j,bi,bj) .GT. 115. .AND.                if ( xC(i,j,bi,bj) .GT. 120. .AND.
547       &             xC(i,j,bi,bj) .LT. 130. .AND.       &             xC(i,j,bi,bj) .LT. 130. .AND.
548       &             yC(i,j,bi,bj) .GT. -10. .AND.       &             yC(i,j,bi,bj) .GT. -10. .AND.
549       &             yC(i,j,bi,bj) .LT.  10. ) then       &             yC(i,j,bi,bj) .LT.  10. ) then
550                   wtp(i,j,bi,bj)  = 0.                   wtp(i,j,bi,bj)  = wtp(i,j,bi,bj)*100.
551                   wers(i,j,bi,bj) = 0.                   wers(i,j,bi,bj) = wers(i,j,bi,bj)*100.
                  wgfo(i,j,bi,bj) = 0.  
552                endif                endif
553  cph-indonesian)  cph-indonesian)
 cph-medit(  
               if ( ( xC(i,j,bi,bj) .GT. 355. .AND.  
      &               xC(i,j,bi,bj) .LT. 360. .AND.  
      &               yC(i,j,bi,bj) .GT.  30. .AND.  
      &               yC(i,j,bi,bj) .LT.  48. )  
      &             .OR.  
      &             ( xC(i,j,bi,bj) .GT.   0. .AND.  
      &               xC(i,j,bi,bj) .LT.  39. .AND.  
      &               yC(i,j,bi,bj) .GT.  30. .AND.  
      &               yC(i,j,bi,bj) .LT.  48. ) ) then  
                   wtp(i,j,bi,bj)  = wtp(i,j,bi,bj) *10.  
                   wers(i,j,bi,bj) = wers(i,j,bi,bj)*10.  
                   wgfo(i,j,bi,bj) = wgfo(i,j,bi,bj)*10.  
               endif  
 cph-medit)  
554              enddo              enddo
555            enddo            enddo
556          enddo          enddo
# Line 623  c--           Test for missing values. Line 576  c--           Test for missing values.
576                if (wscatx(i,j,bi,bj) .lt. -9900.) then                if (wscatx(i,j,bi,bj) .lt. -9900.) then
577                  wscatx(i,j,bi,bj) = 0. _d 0                  wscatx(i,j,bi,bj) = 0. _d 0
578                endif                endif
 c--           Rescale dyn -> N/M^2  
579                wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)                wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)
 c--           Missing values over water should have larger errors  
               if ( wscatx(i,j,bi,bj).EQ.0. .AND.  
      &             maskW(i,j,k,bi,bj).NE.0. )  
      &             wscatx(i,j,bi,bj) = 4.*wtau0  
 c--           Cut off extreme values  
               if ( wscatx(i,j,bi,bj).GT.0.15 )  
      &             wscatx(i,j,bi,bj) = 0.15  
 c--           Set mimimum background  
580                wscatx(i,j,bi,bj) = max(wscatx(i,j,bi,bj),wtau0)                wscatx(i,j,bi,bj) = max(wscatx(i,j,bi,bj),wtau0)
581                wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)*maskW(i,j,k,bi,bj)                wscatx(i,j,bi,bj) = wscatx(i,j,bi,bj)*maskw(i,j,k,bi,bj)
582       &                            *frame(i,j)       &                            *frame(i,j)
 c  
583                if (wscaty(i,j,bi,bj) .lt. -9900.) then                if (wscaty(i,j,bi,bj) .lt. -9900.) then
584                  wscaty(i,j,bi,bj) = 0. _d 0                  wscaty(i,j,bi,bj) = 0. _d 0
585                endif                endif
 c--           Rescale dyn -> N/M^2  
586                wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)                wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)
 c--           Missing values over water should have larger errors  
               if ( wscaty(i,j,bi,bj).EQ.0. .AND.  
      &             maskS(i,j,k,bi,bj).NE.0. )  
      &             wscaty(i,j,bi,bj) = 4.*wtau0  
 c--           Cut off extreme values  
               if ( wscaty(i,j,bi,bj).GT.0.15 )  
      &             wscaty(i,j,bi,bj) = 0.15  
 c--           Set mimimum background  
587                wscaty(i,j,bi,bj) = max(wscaty(i,j,bi,bj),wtau0)                wscaty(i,j,bi,bj) = max(wscaty(i,j,bi,bj),wtau0)
588                wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)*maskS(i,j,k,bi,bj)                wscaty(i,j,bi,bj) = wscaty(i,j,bi,bj)*masks(i,j,k,bi,bj)
589       &                            *frame(i,j)       &                            *frame(i,j)
590              enddo              enddo
591            enddo            enddo
# Line 692  c--           Test for missing values. Line 626  c--           Test for missing values.
626        enddo        enddo
627  #endif  #endif
628    
629  #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION))  #if (defined (ALLOW_USTRESS_COST_CONTRIBUTION) || defined (ALLOW_USTRESS_CONTROL))
630        nnz   =   1        nnz   =   1
631  ce      irec  =   2  ce      irec  =   2
632  ce(   due to Patrick's processing:  ce(   due to Patrick's processing:
# Line 700  ce(   due to Patrick's processing: Line 634  ce(   due to Patrick's processing:
634  ce)  ce)
635        if ( tauu_errfile .NE. ' ' ) then        if ( tauu_errfile .NE. ' ' ) then
636           call mdsreadfield( tauu_errfile, cost_iprec, cost_yftype, nnz,           call mdsreadfield( tauu_errfile, cost_iprec, cost_yftype, nnz,
637       &                   wtauu, irec, mythid )       &        wtauu, irec, mythid )
638        endif        endif
639    
640        do bj = jtlo,jthi        do bj = jtlo,jthi
# Line 712  c--           Test for missing values. Line 646  c--           Test for missing values.
646                if (wtauu(i,j,bi,bj) .lt. -9900.) then                if (wtauu(i,j,bi,bj) .lt. -9900.) then
647                  wtauu(i,j,bi,bj) = 0. _d 0                  wtauu(i,j,bi,bj) = 0. _d 0
648                endif                endif
 c--           Rescale dyn -> N/M^2  
649                wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*0.1                wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*0.1
 c--           Missing values over water should have larger errors  
               if ( wtauu(i,j,bi,bj).EQ.0. .AND.  
      &             maskW(i,j,k,bi,bj).NE.0. )  
      &             wtauu(i,j,bi,bj) = 4.*wtau0  
 c--           Cut off extreme values  
               if ( wtauu(i,j,bi,bj).GT.0.12 )  
      &             wtauu(i,j,bi,bj) = 0.12  
 c--           Set mimimum background  
650                wtauu(i,j,bi,bj) = max(wtauu(i,j,bi,bj),wtau0)                wtauu(i,j,bi,bj) = max(wtauu(i,j,bi,bj),wtau0)
651                wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*maskW(i,j,k,bi,bj)                wtauu(i,j,bi,bj) = wtauu(i,j,bi,bj)*maskw(i,j,k,bi,bj)
652       &                            *frame(i,j)       &                            *frame(i,j)
653  cph(  cph(
654  cph              wtauu2(i,j,bi,bj) = 2.*wtau0*maskW(i,j,k,bi,bj)*frame(i,j)                wtauu2(i,j,bi,bj) = wtau0*maskW(i,j,k,bi,bj)*frame(i,j)
               wtauu2(i,j,bi,bj) = wtauu(i,j,bi,bj)  
655  cph)  cph)
656               enddo               enddo
657            enddo            enddo
658          enddo          enddo
659        enddo        enddo
660    
661  #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION))  #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION) || defined (ALLOW_UWIND_CONTROL))
662    
663        nnz   =   1        nnz   =   1
664  ce      irec  =   2  ce      irec  =   2
# Line 743  ce(   due to Patrick's processing: Line 667  ce(   due to Patrick's processing:
667  ce)  ce)
668        if ( uwind_errfile .NE. ' ' ) then        if ( uwind_errfile .NE. ' ' ) then
669           call mdsreadfield( uwind_errfile, cost_iprec, cost_yftype, nnz,           call mdsreadfield( uwind_errfile, cost_iprec, cost_yftype, nnz,
670       &                   wuwind, irec, mythid )       &        wuwind, irec, mythid )
671        endif        endif
672    
673        do bj = jtlo,jthi        do bj = jtlo,jthi
# Line 757  c--           Test for missing values. Line 681  c--           Test for missing values.
681                endif                endif
682                wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)                wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)
683                wuwind(i,j,bi,bj) = max(wuwind(i,j,bi,bj),wwind0)                wuwind(i,j,bi,bj) = max(wuwind(i,j,bi,bj),wwind0)
684                wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)*maskC(i,j,k,bi,bj)                wuwind(i,j,bi,bj) = wuwind(i,j,bi,bj)*maskc(i,j,k,bi,bj)
685       &                            *frame(i,j)       &                            *frame(i,j)
686               enddo               enddo
687            enddo            enddo
# Line 766  c--           Test for missing values. Line 690  c--           Test for missing values.
690  #endif  #endif
691    
692  c--   Read meridional wind stress variance.  c--   Read meridional wind stress variance.
693  #if (defined (ALLOW_VSTRESS_COST_CONTRIBUTION))  #if (defined (ALLOW_VSTRESS_COST_CONTRIBUTION) || defined (ALLOW_VSTRESS_CONTROL))
694        nnz   =   1        nnz   =   1
695  ce      irec  =   2  ce      irec  =   2
696  ce(   due to Patrick's processing:  ce(   due to Patrick's processing:
# Line 775  ce) Line 699  ce)
699    
700        if ( tauv_errfile .NE. ' ' ) then        if ( tauv_errfile .NE. ' ' ) then
701           call mdsreadfield( tauv_errfile, cost_iprec, cost_yftype, nnz,           call mdsreadfield( tauv_errfile, cost_iprec, cost_yftype, nnz,
702       &                   wtauv, irec, mythid )       &        wtauv, irec, mythid )
703        endif        endif
704    
705        do bj = jtlo,jthi        do bj = jtlo,jthi
# Line 786  c--           Test for missing values. Line 710  c--           Test for missing values.
710                if (wtauv(i,j,bi,bj) .lt. -9900.) then                if (wtauv(i,j,bi,bj) .lt. -9900.) then
711                  wtauv(i,j,bi,bj) = 0. _d 0                  wtauv(i,j,bi,bj) = 0. _d 0
712                endif                endif
 c--           Rescape dyn -> dyn  
713                wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*0.1                wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*0.1
 c--           Missing values over water should have larger errors  
               if ( wtauv(i,j,bi,bj).EQ.0. .AND.  
      &             maskS(i,j,k,bi,bj).NE.0. )  
      &             wtauv(i,j,bi,bj) = 4.*wtau0  
 c--           Cut off extreme values  
               if ( wtauv(i,j,bi,bj).GT.0.12 )  
      &             wtauv(i,j,bi,bj) = 0.12  
 c--           Set mimimum background  
714                wtauv(i,j,bi,bj) = max(wtauv(i,j,bi,bj),wtau0)                wtauv(i,j,bi,bj) = max(wtauv(i,j,bi,bj),wtau0)
715                wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*maskS(i,j,k,bi,bj)                wtauv(i,j,bi,bj) = wtauv(i,j,bi,bj)*masks(i,j,k,bi,bj)
716       &                            *frame(i,j)       &                            *frame(i,j)
717  cph(  cph(
718  cph              wtauv2(i,j,bi,bj) = 2.*wtau0*maskS(i,j,k,bi,bj)*frame(i,j)                wtauv2(i,j,bi,bj) = wtau0*maskS(i,j,k,bi,bj)*frame(i,j)
               wtauv2(i,j,bi,bj) = wtauv(i,j,bi,bj)  
719  cph)  cph)
720              enddo              enddo
721            enddo            enddo
722          enddo          enddo
723        enddo        enddo
724    
725  #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION))  #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION) || defined (ALLOW_VWIND_CONTROL))
726    
727        nnz   =   1        nnz   =   1
728  ce      irec  =   2  ce      irec  =   2
# Line 818  ce) Line 732  ce)
732    
733        if ( vwind_errfile .NE. ' ' ) then        if ( vwind_errfile .NE. ' ' ) then
734           call mdsreadfield( vwind_errfile, cost_iprec, cost_yftype, nnz,           call mdsreadfield( vwind_errfile, cost_iprec, cost_yftype, nnz,
735       &                   wvwind, irec, mythid )       &        wvwind, irec, mythid )
736        endif        endif
737    
738        do bj = jtlo,jthi        do bj = jtlo,jthi
# Line 831  c--           Test for missing values. Line 745  c--           Test for missing values.
745                endif                endif
746                wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)                wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)
747                wvwind(i,j,bi,bj) = max(wvwind(i,j,bi,bj),wwind0)                wvwind(i,j,bi,bj) = max(wvwind(i,j,bi,bj),wwind0)
748                wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)*maskC(i,j,k,bi,bj)                wvwind(i,j,bi,bj) = wvwind(i,j,bi,bj)*maskc(i,j,k,bi,bj)
749       &                            *frame(i,j)       &                            *frame(i,j)
750               enddo              enddo
751            enddo            enddo
752          enddo          enddo
753        enddo        enddo
754  #endif  #endif
755    
756  #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION))  #if (defined (ALLOW_HFLUX_COST_CONTRIBUTION) || defined (ALLOW_HFLUX_CONTROL))
757  c--   Read heat flux flux variance.  c--   Read heat flux flux variance.
758        nnz   =  1        nnz   =  1
759  c--   First  record in data file:  mean field.  c--   First  record in data file:  mean field.
# Line 849  ce(   due to Patrick's processing: Line 763  ce(   due to Patrick's processing:
763        irec  = 1        irec  = 1
764  ce)  ce)
765        if ( hflux_errfile .NE. ' ' ) then        if ( hflux_errfile .NE. ' ' ) then
766           call mdsreadfield( hflux_errfile, cost_iprec, cost_yftype, nnz,           call mdsreadfield( hflux_errfile, cost_iprec, cost_yftype, nnz,
767       &                   whflux, irec, mythid )       &        whflux, irec, mythid )
768        endif        endif
769    
770        do bj = jtlo,jthi        do bj = jtlo,jthi
# Line 868  c--           Data are in units of W/m** Line 782  c--           Data are in units of W/m**
782                whfluxm(i,j,bi,bj) = max(whfluxm(i,j,bi,bj),whflux0m)                whfluxm(i,j,bi,bj) = max(whfluxm(i,j,bi,bj),whflux0m)
783       &                            *frame(i,j)       &                            *frame(i,j)
784  cph(  cph(
785  cph              whflux2(i,j,bi,bj) = 2.*whflux0*frame(i,j)                whflux2(i,j,bi,bj) = whflux0*frame(i,j)
               whflux2(i,j,bi,bj) = whflux(i,j,bi,bj)  
786  cph)  cph)
787              enddo              enddo
788            enddo            enddo
789          enddo          enddo
790        enddo        enddo
791  #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION))  #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION) || defined (ALLOW_ATEMP_CONTROL))
792  c--   Read atmos. temp. variance.  c--   Read atmos. temp. variance.
793        nnz   =  1        nnz   =  1
794  c--   First  record in data file:  mean field.  c--   First  record in data file:  mean field.
# Line 886  ce(   due to Patrick's processing: Line 799  ce(   due to Patrick's processing:
799  ce)  ce)
800        if ( atemp_errfile .NE. ' ' ) then        if ( atemp_errfile .NE. ' ' ) then
801           call mdsreadfield( atemp_errfile, cost_iprec, cost_yftype, nnz,           call mdsreadfield( atemp_errfile, cost_iprec, cost_yftype, nnz,
802       &                   watemp, irec, mythid )       &        watemp, irec, mythid )
803        endif        endif
804    
805        do bj = jtlo,jthi        do bj = jtlo,jthi
# Line 897  c--           Test for missing values. Line 810  c--           Test for missing values.
810                if (watemp(i,j,bi,bj) .lt. -9900.) then                if (watemp(i,j,bi,bj) .lt. -9900.) then
811                  watemp(i,j,bi,bj) = 0. _d 0                  watemp(i,j,bi,bj) = 0. _d 0
812                endif                endif
813  c--           Data are in units of deg.  c--           Data are in units of W/m**2.
814                watemp(i,j,bi,bj) = watemp(i,j,bi,bj)                watemp(i,j,bi,bj) = watemp(i,j,bi,bj)
815                watemp(i,j,bi,bj) = max(watemp(i,j,bi,bj),watemp0)                watemp(i,j,bi,bj) = max(watemp(i,j,bi,bj),watemp0)
816       &                            *frame(i,j)       &                            *frame(i,j)
# Line 905  c--           Data are in units of deg. Line 818  c--           Data are in units of deg.
818            enddo            enddo
819          enddo          enddo
820        enddo        enddo
821    
822    
823  #endif  #endif
824    
825  #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION))  #if (defined (ALLOW_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SFLUX_CONTROL))
826  c--   Read salt flux variance. Second read: data in units of m/s.  c--   Read salt flux variance. Second read: data in units of m/s.
827        nnz   =  1        nnz   =  1
828  c--   First  record in data file:  mean field.  c--   First  record in data file:  mean field.
# Line 917  ce(   due to Patrick's processing: Line 832  ce(   due to Patrick's processing:
832        irec  = 1        irec  = 1
833  ce)  ce)
834        if ( sflux_errfile .NE. ' ' ) then        if ( sflux_errfile .NE. ' ' ) then
835           call mdsreadfield( sflux_errfile, cost_iprec, cost_yftype, nnz,           call mdsreadfield( sflux_errfile, cost_iprec, cost_yftype, nnz,
836       &                   wsflux, irec, mythid )       &        wsflux, irec, mythid )
837        endif        endif
838    
839        do bj = jtlo,jthi        do bj = jtlo,jthi
# Line 936  c--           Data are in units of m/s. Line 851  c--           Data are in units of m/s.
851                wsfluxm(i,j,bi,bj) = max(wsfluxm(i,j,bi,bj),wsflux0m)                wsfluxm(i,j,bi,bj) = max(wsfluxm(i,j,bi,bj),wsflux0m)
852       &                            *frame(i,j)       &                            *frame(i,j)
853  cph(  cph(
854  cph              wsflux2(i,j,bi,bj) = 2.*wsflux0*frame(i,j)                wsflux2(i,j,bi,bj) = wsflux0*frame(i,j)
               wsflux2(i,j,bi,bj) = wsflux(i,j,bi,bj)  
855  cph)  cph)
856              enddo              enddo
857            enddo            enddo
858          enddo          enddo
859        enddo        enddo
860  #elif (defined (ALLOW_AQH_COST_CONTRIBUTION))  #elif (defined (ALLOW_AQH_COST_CONTRIBUTION) || defined (ALLOW_AQH_CONTROL))
861  c--   Secific humid. variance. Second read: data in units of m/s.  c--   Secific humid. variance. Second read: data in units of m/s.
862        nnz   =  1        nnz   =  1
863  c--   First  record in data file:  mean field.  c--   First  record in data file:  mean field.
# Line 953  ce(   due to Patrick's processing: Line 867  ce(   due to Patrick's processing:
867        irec  = 1        irec  = 1
868  ce)  ce)
869        if ( aqh_errfile .NE. ' ' ) then        if ( aqh_errfile .NE. ' ' ) then
870           call mdsreadfield( aqh_errfile, cost_iprec, cost_yftype, nnz,           call mdsreadfield( aqh_errfile, cost_iprec, cost_yftype, nnz,
871       &                   waqh, irec, mythid )       &        waqh, irec, mythid )
872        endif        endif
873    
874        do bj = jtlo,jthi        do bj = jtlo,jthi
# Line 965  c--           Test for missing values. Line 879  c--           Test for missing values.
879                if (waqh(i,j,bi,bj) .lt. -9900.) then                if (waqh(i,j,bi,bj) .lt. -9900.) then
880                  waqh(i,j,bi,bj) = 0. _d 0                  waqh(i,j,bi,bj) = 0. _d 0
881                endif                endif
882  c--           Data are in units of  c--           Data are in units of m/s.
883                waqh(i,j,bi,bj) = waqh(i,j,bi,bj)                waqh(i,j,bi,bj) = waqh(i,j,bi,bj)
884                waqh(i,j,bi,bj) = max(waqh(i,j,bi,bj),waqh0)                waqh(i,j,bi,bj) = max(waqh(i,j,bi,bj),waqh0)
885       &                            *frame(i,j)       &                            *frame(i,j)
# Line 975  c--           Data are in units of Line 889  c--           Data are in units of
889        enddo        enddo
890  #endif  #endif
891    
892    #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  c--   Units have to be checked!  c--   Units have to be checked!
989        do bj = jtlo,jthi        do bj = jtlo,jthi
990          do bi = itlo,ithi          do bi = itlo,ithi
# Line 986  c--   Units have to be checked! Line 996  c--   Units have to be checked!
996                if (wers(i,j,bi,bj) .ne. 0.) then                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)                  wers(i,j,bi,bj) = 1./wers(i,j,bi,bj)/wers(i,j,bi,bj)
998                endif                endif
               if (wgfo(i,j,bi,bj) .ne. 0.) then  
                 wgfo(i,j,bi,bj) = 1./wgfo(i,j,bi,bj)/wgfo(i,j,bi,bj)  
               endif  
 cph(  
 cph sst, sss: reduce weights by factor 2  
999                if (wsst(i,j,bi,bj) .ne. 0.) then                if (wsst(i,j,bi,bj) .ne. 0.) then
1000                  wsst(i,j,bi,bj) = 1./wsst(i,j,bi,bj)/wsst(i,j,bi,bj)/2.                  wsst(i,j,bi,bj) = 1./wsst(i,j,bi,bj)/wsst(i,j,bi,bj)
1001                endif                endif
1002                if (wsss(i,j,bi,bj) .ne. 0.) then                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)/2.                  wsss(i,j,bi,bj) = 1./wsss(i,j,bi,bj)/wsss(i,j,bi,bj)
1004                endif                endif
 cph)  
1005                if (wp(i,j,bi,bj) .ne. 0.) then                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)                  wp(i,j,bi,bj) = 1./wp(i,j,bi,bj)/wp(i,j,bi,bj)
1007                endif                endif
1008                if (wtauu(i,j,bi,bj) .ne. 0.) then                if (wtauu(i,j,bi,bj) .ne. 0.) then
1009                  wtauu(i,j,bi,bj) =                  wtauu(i,j,bi,bj) = 1./wtauu(i,j,bi,bj)/wtauu(i,j,bi,bj)
      &               1./wtauu(i,j,bi,bj)/wtauu(i,j,bi,bj)  
1010                else                else
1011                  wtauu(i,j,bi,bj) = 0.0 _d 0                  wtauu(i,j,bi,bj) = 0.0 _d 0
1012                endif                endif
# Line 1020  cph) Line 1023  cph)
1023                  wscatx(i,j,bi,bj) = 0.0 _d 0                  wscatx(i,j,bi,bj) = 0.0 _d 0
1024                endif                endif
1025                if (wtauv(i,j,bi,bj) .ne. 0.) then                if (wtauv(i,j,bi,bj) .ne. 0.) then
1026                  wtauv(i,j,bi,bj) =                  wtauv(i,j,bi,bj) = 1./wtauv(i,j,bi,bj)/wtauv(i,j,bi,bj)
      &               1./wtauv(i,j,bi,bj)/wtauv(i,j,bi,bj)  
1027                else                else
1028                  wtauv(i,j,bi,bj) = 0.0 _d 0                  wtauv(i,j,bi,bj) = 0.0 _d 0
1029                endif                endif
# Line 1061  cph) Line 1063  cph)
1063                else                else
1064                  wsfluxm(i,j,bi,bj) = 0.0 _d 0                  wsfluxm(i,j,bi,bj) = 0.0 _d 0
1065                endif                endif
 cph)  
1066                if (wuwind(i,j,bi,bj) .ne. 0.) then                if (wuwind(i,j,bi,bj) .ne. 0.) then
1067                  wuwind(i,j,bi,bj) =                  wuwind(i,j,bi,bj) =
1068       &                1./wuwind(i,j,bi,bj)/wuwind(i,j,bi,bj)       &                1./wuwind(i,j,bi,bj)/wuwind(i,j,bi,bj)
# Line 1086  cph) Line 1087  cph)
1087                else                else
1088                  waqh(i,j,bi,bj) = 0.0 _d 0                  waqh(i,j,bi,bj) = 0.0 _d 0
1089                endif                endif
1090                  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    
1109                if (wsfluxmm(bi,bj).ne.0.)                if (wsfluxmm(bi,bj).ne.0.)
1110       &             wsfluxmm(bi,bj) = 1./wsfluxmm(bi,bj)*wsfluxmm(bi,bj)       &             wsfluxmm(bi,bj) = 1./wsfluxmm(bi,bj)*wsfluxmm(bi,bj)
1111                if (whfluxmm(bi,bj).ne.0.)                if (whfluxmm(bi,bj).ne.0.)
1112       &             whfluxmm(bi,bj) = 1./whfluxmm(bi,bj)*whfluxmm(bi,bj)       &             whfluxmm(bi,bj) = 1./whfluxmm(bi,bj)*whfluxmm(bi,bj)
1113    
1114  cph(  cph(
1115                if (whflux2(i,j,bi,bj) .ne. 0.) then                if (whflux2(i,j,bi,bj) .ne. 0.) then
1116                   whflux2(i,j,bi,bj) =                   whflux2(i,j,bi,bj) =
# Line 1105  cph( Line 1125  cph(
1125                   wsflux2(i,j,bi,bj) = 0.0 _d 0                   wsflux2(i,j,bi,bj) = 0.0 _d 0
1126                endif                endif
1127                if (wtauu2(i,j,bi,bj) .ne. 0.) then                if (wtauu2(i,j,bi,bj) .ne. 0.) then
1128                   wtauu2(i,j,bi,bj) =                   wtauu2(i,j,bi,bj) =
1129       &                1./wtauu2(i,j,bi,bj)/wtauu2(i,j,bi,bj)       &                1./wtauu2(i,j,bi,bj)/wtauu2(i,j,bi,bj)
1130                else                else
1131                   wtauu2(i,j,bi,bj) = 0.0 _d 0                   wtauu2(i,j,bi,bj) = 0.0 _d 0
# Line 1120  cph) Line 1140  cph)
1140              enddo              enddo
1141            enddo            enddo
1142    
 cph(  
 cph reduce wtheta, wsalt by factor 2.  
1143            do k = 1,nr            do k = 1,nr
1144              if (wtheta(k,bi,bj) .ne. 0.) then              if (wtheta(k,bi,bj) .ne. 0.) then
1145                wtheta(k,bi,bj) = ratio/wtheta(k,bi,bj)/wtheta(k,bi,bj)/2.                wtheta(k,bi,bj) = ratio/wtheta(k,bi,bj)/wtheta(k,bi,bj)
1146              else              else
1147                wtheta(k,bi,bj) = 0.0 _d 0                wtheta(k,bi,bj) = 0.0 _d 0
1148              endif              endif
1149              if (wsalt(k,bi,bj) .ne. 0.) then              if (wsalt(k,bi,bj) .ne. 0.) then
1150                wsalt(k,bi,bj) = ratio/wsalt(k,bi,bj)/wsalt(k,bi,bj)/2.                wsalt(k,bi,bj) = ratio/wsalt(k,bi,bj)/wsalt(k,bi,bj)
1151              else              else
1152                wsalt(k,bi,bj) = 0.0 _d 0                wsalt(k,bi,bj) = 0.0 _d 0
1153              endif              endif
1154            enddo            enddo
 cph)  
1155    
1156  #ifdef ALLOW_OBCS_COST_CONTRIBUTION  #ifdef ALLOW_OBCS_COST_CONTRIBUTION
1157            do iobcs = 1,nobcs            do iobcs = 1,nobcs
# Line 1178  cph) Line 1195  cph)
1195          enddo          enddo
1196        enddo        enddo
1197    
1198  #if   (defined (ALLOW_HFLUX_COST_CONTRIBUTION))        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        call active_write_xy_loc( 'whflux', whflux, 1, 0, mythid, dummy)        call active_write_xy_loc( 'whflux', whflux, 1, 0, mythid, dummy)
1271        call active_write_xy_loc( 'whflux2', whflux2, 1, 0, mythid, dummy)        call active_write_xy_loc( 'whflux2', whflux2, 1, 0, mythid, dummy)
1272  #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION))  #elif (defined (ALLOW_ATEMP_COST_CONTRIBUTION) || defined (ALLOW_ATEMP_CONTROL))
1273        call active_write_xy_loc( 'watemp', watemp, 1, 0, mythid, dummy)        call active_write_xy_loc( 'watemp', watemp, 1, 0, mythid, dummy)
1274  #endif  #endif
1275    
1276  #if   (defined (ALLOW_SFLUX_COST_CONTRIBUTION))  #if   (defined (ALLOW_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SFLUX_CONTROL))
1277        call active_write_xy_loc( 'wsflux', wsflux, 1, 0, mythid, dummy)        call active_write_xy_loc( 'wsflux', wsflux, 1, 0, mythid, dummy)
1278        call active_write_xy_loc( 'wsflux2', wsflux2, 1, 0, mythid, dummy)        call active_write_xy_loc( 'wsflux2', wsflux2, 1, 0, mythid, dummy)
1279  #elif (defined (ALLOW_AQH_COST_CONTRIBUTION))  #elif (defined (ALLOW_AQH_COST_CONTRIBUTION) || defined (ALLOW_AQH_CONTROL))
1280        call active_write_xy_loc( 'waqh', waqh, 1, 0, mythid, dummy)        call active_write_xy_loc( 'waqh', waqh, 1, 0, mythid, dummy)
1281  #endif  #endif
1282    
1283  #if   (defined (ALLOW_USTRESS_COST_CONTRIBUTION))  #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        call active_write_xy_loc( 'wtauu', wtauu,   1, 0, mythid, dummy)        call active_write_xy_loc( 'wtauu', wtauu,   1, 0, mythid, dummy)
1297        call active_write_xy_loc( 'wtauu2', wtauu2,   1, 0, mythid, dummy)        call active_write_xy_loc( 'wtauu2', wtauu2,   1, 0, mythid, dummy)
1298  #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION))  #elif (defined (ALLOW_UWIND_COST_CONTRIBUTION) || defined (ALLOW_UWIND_CONTROL))
1299        call active_write_xy_loc( 'wuwind', wuwind, 1, 0, mythid, dummy)        call active_write_xy_loc( 'wuwind', wuwind, 1, 0, mythid, dummy)
1300  #endif  #endif
1301    
1302  #if   (defined (ALLOW_VSTRESS_COST_CONTRIBUTION))  #if   (defined (ALLOW_VSTRESS_COST_CONTRIBUTION) || defined (ALLOW_VSTRESS_CONTROL))
1303        call active_write_xy_loc( 'wtauv', wtauv,   1, 0, mythid, dummy)        call active_write_xy_loc( 'wtauv', wtauv,   1, 0, mythid, dummy)
1304        call active_write_xy_loc( 'wtauv2', wtauv2,   1, 0, mythid, dummy)        call active_write_xy_loc( 'wtauv2', wtauv2,   1, 0, mythid, dummy)
1305  #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION))  #elif (defined (ALLOW_VWIND_COST_CONTRIBUTION) || defined (ALLOW_VWIND_CONTROL))
1306        call active_write_xy_loc( 'wvwind', wvwind, 1, 0, mythid, dummy)        call active_write_xy_loc( 'wvwind', wvwind, 1, 0, mythid, dummy)
1307  #endif  #endif
1308    
1309    #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  #ifdef ALLOW_OBCSN_COST_CONTRIBUTION  #ifdef ALLOW_OBCSN_COST_CONTRIBUTION
1318  #endif  #endif
1319    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22