tutorial_global_oce_optim/code_ad/cost_weights.F

C $Header: /u/gcmpack/MITgcm/verification/tutorial_global_oce_optim/code_ad/cost_weights.F,v 1.2 2008/06/23 18:38:34 dfer Exp $
C $Name:  $

#include "COST_CPPOPTIONS.h"

      SUBROUTINE COST_WEIGHTS( myThid )

C     ==================================================================
C     SUBROUTINE COST_WEIGHTS
C     ==================================================================
C
C     o Set weights used in the cost function and in the
C       normalization of the sensitivities when ALLOW_NON_DIMENSIONAL

      IMPLICIT NONE

C     == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"

#include "ctrl.h"
#include "ctrl_weights.h"
#include "cost.h"

C     == routine arguments ==
      INTEGER  myThid

C     == Functions ==
      INTEGER  MDS_RECLEN
      EXTERNAL MDS_RECLEN

C     == local variables ==
      INTEGER bi,bj
      INTEGER i,j,k
      INTEGER itlo,ithi,jtlo,jthi
      INTEGER jMin,jMax,iMin,iMax
      INTEGER iUnit, length_of_rec

      _RL dummy
      _RL wti(Nr)
      REAL*8 tmpwti(Nr)
      CHARACTER*(MAX_LEN_MBUF) msgBuf

C     == end of interface ==

      jtlo = myByLo(myThid)
      jthi = myByHi(myThid)
      itlo = myBxLo(myThid)
      ithi = myBxHi(myThid)
      iMin = 1-OLx
      iMax = sNx+OLx
      jMin = 1-OLy
      jMax = sNy+OLy

C--   Initialize variance (weight) fields.
      DO k = 1,Nr
         wti(k) = 0. _d 0
      ENDDO
      DO bj = jtlo,jthi
        DO bi = itlo,ithi
          DO j = jMin,jMax
            DO i = iMin,iMax
              whfluxm(i,j,bi,bj)= 0. _d 0
            ENDDO
          ENDDO
          DO k = 1,Nr
             wunit(k,bi,bj)  = 1. _d 0
             wtheta(k,bi,bj) = 0. _d 0
             wsalt(k,bi,bj)  = 0. _d 0
          ENDDO
        ENDDO
      ENDDO

C--   Read error information and set up weight matrices.

#ifdef ALLOW_COST_TEMP
C  Temperature weights for cost function
       _BEGIN_MASTER(myThid)
       CALL MDSFINDUNIT( iUnit, myThid )
       length_of_rec = MDS_RECLEN( precFloat64, Nr, myThid )
       OPEN( iUnit, FILE='Err_levitus_15layer.bin', STATUS='OLD',
     &       FORM='UNFORMATTED',ACCESS='DIRECT',RECL=length_of_rec )
       READ(iUnit,rec=1) tmpwti
       CLOSE(iUnit)
#ifdef _BYTESWAPIO
       CALL MDS_BYTESWAPR8( Nr, tmpwti )
#endif
       _END_MASTER(myThid)
       _BARRIER

       DO k=1,Nr
         wti(k) = tmpwti(k)
       ENDDO
       WRITE(msgBuf,'(3A)') 'S/R COST_WEIGHTS:',
     &  ' Temperature weights loaded from: ','Err_levitus_15layer.bin'
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                     SQUEEZE_RIGHT , myThid )

      print*,'Weights for temperature: wti', (wti(k),k=1,nr)

      DO bj = jtlo,jthi
        DO bi = itlo,ithi
          DO k = 1, Nr
               wtheta(k,bi,bj) = 1. _d 0/wti(k)/wti(k)
          ENDDO
        ENDDO
      ENDDO
#endif

C--   Then the hflux weights :

#if (defined (ALLOW_COST_HFLUXM) || defined (ALLOW_HFLUXM_CONTROL))
      CALL READ_REC_3D_RL( 'Err_hflux.bin', precFloat64, 1,
     &                      whfluxm, 1, 0, myThid )
      _EXCH_XY_RL(whfluxm   , myThid )
      DO bj = jtlo,jthi
        DO bi = itlo,ithi
          DO j = jMin,jMax
            DO i = iMin,iMax
             print*,'Uncertainties for Heat Flux',i,j,whfluxm(i,j,bi,bj)
             IF (whfluxm(i,j,bi,bj) .NE. 0. _d 0) THEN
                 whfluxm(i,j,bi,bj) = 1. _d 0 /whfluxm(i,j,bi,bj)
     &                                        /whfluxm(i,j,bi,bj)
             ELSE
                 whfluxm(i,j,bi,bj) = 1. _d 0
             ENDIF
            ENDDO
          ENDDO
        ENDDO
      ENDDO
#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
      CALL ACTIVE_WRITE_XY('whfluxm',whfluxm,1,0,myThid,dummy)
#endif
#endif
      RETURN
      END
1	jmc	1.3	C $Header: /u/gcmpack/MITgcm/verification/tutorial_global_oce_optim/code_ad/cost_weights.F,v 1.2 2008/06/23 18:38:34 dfer Exp $
2	dfer	1.1	C $Name: $
3
4			#include "COST_CPPOPTIONS.h"
5
6	dfer	1.2	SUBROUTINE COST_WEIGHTS( myThid )
7	dfer	1.1
8	jmc	1.3	C ==================================================================
9			C SUBROUTINE COST_WEIGHTS
10			C ==================================================================
11			C
12			C o Set weights used in the cost function and in the
13			C normalization of the sensitivities when ALLOW_NON_DIMENSIONAL
14	dfer	1.1
15	dfer	1.2	IMPLICIT NONE
16	dfer	1.1
17	jmc	1.3	C == global variables ==
18	dfer	1.1	#include "EEPARAMS.h"
19			#include "SIZE.h"
20			#include "PARAMS.h"
21			#include "GRID.h"
22
23			#include "ctrl.h"
24			#include "ctrl_weights.h"
25			#include "cost.h"
26
27	jmc	1.3	C == routine arguments ==
28	dfer	1.2	INTEGER myThid
29	dfer	1.1
30	jmc	1.3	C == Functions ==
31			INTEGER MDS_RECLEN
32			EXTERNAL MDS_RECLEN
33	dfer	1.1
34	jmc	1.3	C == local variables ==
35	dfer	1.2	INTEGER bi,bj
36			INTEGER i,j,k
37			INTEGER itlo,ithi,jtlo,jthi
38			INTEGER jMin,jMax,iMin,iMax
39	jmc	1.3	INTEGER iUnit, length_of_rec
40	dfer	1.1
41			_RL dummy
42	dfer	1.2	_RL wti(Nr)
43			REAL*8 tmpwti(Nr)
44			CHARACTER*(MAX_LEN_MBUF) msgBuf
45	dfer	1.1
46	jmc	1.3	C == end of interface ==
47	dfer	1.1
48	dfer	1.2	jtlo = myByLo(myThid)
49			jthi = myByHi(myThid)
50			itlo = myBxLo(myThid)
51			ithi = myBxHi(myThid)
52			iMin = 1-OLx
53			iMax = sNx+OLx
54			jMin = 1-OLy
55			jMax = sNy+OLy
56	dfer	1.1
57	jmc	1.3	C-- Initialize variance (weight) fields.
58	dfer	1.2	DO k = 1,Nr
59			wti(k) = 0. _d 0
60			ENDDO
61	dfer	1.1	DO bj = jtlo,jthi
62			DO bi = itlo,ithi
63	dfer	1.2	DO j = jMin,jMax
64			DO i = iMin,iMax
65	dfer	1.1	whfluxm(i,j,bi,bj)= 0. _d 0
66			ENDDO
67			ENDDO
68	dfer	1.2	DO k = 1,Nr
69	dfer	1.1	wunit(k,bi,bj) = 1. _d 0
70			wtheta(k,bi,bj) = 0. _d 0
71			wsalt(k,bi,bj) = 0. _d 0
72			ENDDO
73			ENDDO
74			ENDDO
75
76	jmc	1.3	C-- Read error information and set up weight matrices.
77	dfer	1.1
78			#ifdef ALLOW_COST_TEMP
79	jmc	1.3	C Temperature weights for cost function
80	dfer	1.2	_BEGIN_MASTER(myThid)
81			CALL MDSFINDUNIT( iUnit, myThid )
82	jmc	1.3	length_of_rec = MDS_RECLEN( precFloat64, Nr, myThid )
83			OPEN( iUnit, FILE='Err_levitus_15layer.bin', STATUS='OLD',
84			& FORM='UNFORMATTED',ACCESS='DIRECT',RECL=length_of_rec )
85			READ(iUnit,rec=1) tmpwti
86	dfer	1.2	CLOSE(iUnit)
87			#ifdef _BYTESWAPIO
88			CALL MDS_BYTESWAPR8( Nr, tmpwti )
89			#endif
90			_END_MASTER(myThid)
91			_BARRIER
92
93			DO k=1,Nr
94			wti(k) = tmpwti(k)
95			ENDDO
96			WRITE(msgBuf,'(3A)') 'S/R COST_WEIGHTS:',
97			& ' Temperature weights loaded from: ','Err_levitus_15layer.bin'
98			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
99			& SQUEEZE_RIGHT , myThid )
100
101	dfer	1.1	print*,'Weights for temperature: wti', (wti(k),k=1,nr)
102
103			DO bj = jtlo,jthi
104			DO bi = itlo,ithi
105	dfer	1.2	DO k = 1, Nr
106	dfer	1.1	wtheta(k,bi,bj) = 1. _d 0/wti(k)/wti(k)
107			ENDDO
108			ENDDO
109			ENDDO
110			#endif
111	jmc	1.3
112			C-- Then the hflux weights :
113
114	dfer	1.1	#if (defined (ALLOW_COST_HFLUXM) \|\| defined (ALLOW_HFLUXM_CONTROL))
115	jmc	1.3	CALL READ_REC_3D_RL( 'Err_hflux.bin', precFloat64, 1,
116			& whfluxm, 1, 0, myThid )
117	dfer	1.1	_EXCH_XY_RL(whfluxm , myThid )
118			DO bj = jtlo,jthi
119			DO bi = itlo,ithi
120	dfer	1.2	DO j = jMin,jMax
121			DO i = iMin,iMax
122	dfer	1.1	print*,'Uncertainties for Heat Flux',i,j,whfluxm(i,j,bi,bj)
123			IF (whfluxm(i,j,bi,bj) .NE. 0. _d 0) THEN
124			whfluxm(i,j,bi,bj) = 1. _d 0 /whfluxm(i,j,bi,bj)
125			& /whfluxm(i,j,bi,bj)
126			ELSE
127			whfluxm(i,j,bi,bj) = 1. _d 0
128			ENDIF
129			ENDDO
130			ENDDO
131			ENDDO
132			ENDDO
133			#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
134	jmc	1.3	CALL ACTIVE_WRITE_XY('whfluxm',whfluxm,1,0,myThid,dummy)
135	dfer	1.1	#endif
136			#endif
137	jmc	1.3	RETURN
138			END