tutorial_global_oce_optim/code_oad/cost_weights.F

0e09621e3e Patr*0001 #include "COST_OPTIONS.h"
987a76ae74 Jean*0002 #ifdef ALLOW_CTRL
                0003 # include "CTRL_OPTIONS.h"
                0004 #endif
0e09621e3e Patr*0005
cf705a6c8e Mart*0006 CBOP
                0007 C     !ROUTINE: COST_WEIGHTS
                0008 C     !INTERFACE:
0e09621e3e Patr*0009       SUBROUTINE COST_WEIGHTS( myThid )
cf705a6c8e Mart*0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
                0012 C     | SUBROUTINE COST_WEIGHTS
                0013 C     | o Set weights used in the cost function
                0014 C     *==========================================================*
0e09621e3e Patr*0015 C
cf705a6c8e Mart*0016 C     \ev
                0017 C     !USES:
0e09621e3e Patr*0018       IMPLICIT NONE
cf705a6c8e Mart*0019 C     == Global variables ===
0e09621e3e Patr*0020 #include "EEPARAMS.h"
                0021 #include "SIZE.h"
                0022 #include "PARAMS.h"
                0023 #include "GRID.h"
cf705a6c8e Mart*0024 #ifdef ALLOW_OPENAD
                0025 # include "FFIELDS.h"
                0026 #else
                0027 # include "cost_local.h"
                0028 #endif
0e09621e3e Patr*0029
cf705a6c8e Mart*0030 C     !INPUT/OUTPUT PARAMETERS:
0e09621e3e Patr*0031       INTEGER  myThid
                0032
cf705a6c8e Mart*0033 C     !FUNCTIONS:
0e09621e3e Patr*0034       INTEGER  MDS_RECLEN
                0035       EXTERNAL MDS_RECLEN
                0036
cf705a6c8e Mart*0037 C     !LOCAL VARIABLES:
0e09621e3e Patr*0038       INTEGER bi,bj
                0039       INTEGER i,j,k
                0040       INTEGER iUnit, length_of_rec
                0041
                0042       _RL dummy
                0043       _RL wti(Nr)
                0044       REAL*8 tmpwti(Nr)
                0045       CHARACTER*(MAX_LEN_MBUF) msgBuf
cf705a6c8e Mart*0046 CEOP
0e09621e3e Patr*0047
                0048 C--   Initialize variance (weight) fields.
                0049       DO k = 1,Nr
cf705a6c8e Mart*0050        wti(k) = 0. _d 0
0e09621e3e Patr*0051       ENDDO
cf705a6c8e Mart*0052       DO bj = myByLo(myThid),myByHi(myThid)
                0053        DO bi = myBxLo(myThid),myBxHi(myThid)
                0054         DO j = 1-OLy,sNy+OLy
                0055          DO i = 1-OLx,sNx+OLx
                0056           whfluxm(i,j,bi,bj)= 0. _d 0
                0057          ENDDO
                0058         ENDDO
                0059         DO k = 1,Nr
                0060          wtheta(k,bi,bj) = 0. _d 0
0e09621e3e Patr*0061         ENDDO
cf705a6c8e Mart*0062        ENDDO
0e09621e3e Patr*0063       ENDDO
                0064
                0065 C--   Read error information and set up weight matrices.
                0066
                0067 #ifdef ALLOW_COST_TEMP
                0068 C  Temperature weights for cost function
                0069        _BEGIN_MASTER(myThid)
                0070        CALL MDSFINDUNIT( iUnit, myThid )
                0071        length_of_rec = MDS_RECLEN( precFloat64, Nr, myThid )
                0072        OPEN( iUnit, FILE='Err_levitus_15layer.bin', STATUS='OLD',
                0073      &       FORM='UNFORMATTED',ACCESS='DIRECT',RECL=length_of_rec )
                0074        READ(iUnit,rec=1) tmpwti
                0075        CLOSE(iUnit)
                0076 #ifdef _BYTESWAPIO
                0077        CALL MDS_BYTESWAPR8( Nr, tmpwti )
                0078 #endif
                0079        _END_MASTER(myThid)
                0080        _BARRIER
                0081
                0082        DO k=1,Nr
                0083          wti(k) = tmpwti(k)
                0084        ENDDO
                0085        WRITE(msgBuf,'(3A)') 'S/R COST_WEIGHTS:',
                0086      &  ' Temperature weights loaded from: ','Err_levitus_15layer.bin'
                0087        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0088      &                     SQUEEZE_RIGHT , myThid )
                0089
                0090 c     print*,'Weights for temperature: wti', (wti(k),k=1,nr)
                0091
cf705a6c8e Mart*0092       DO bj = myByLo(myThid),myByHi(myThid)
                0093         DO bi = myBxLo(myThid),myBxHi(myThid)
0e09621e3e Patr*0094           DO k = 1, Nr
                0095                wtheta(k,bi,bj) = 1. _d 0/wti(k)/wti(k)
                0096           ENDDO
                0097         ENDDO
                0098       ENDDO
987a76ae74 Jean*0099 #endif /* ALLOW_COST_TEMP */
0e09621e3e Patr*0100
                0101 C--   Then the hflux weights :
                0102
2e7aec9951 dngo*0103 #ifdef ALLOW_COST_HFLUXM
0e09621e3e Patr*0104       CALL READ_REC_3D_RL( 'Err_hflux.bin', precFloat64, 1,
                0105      &                      whfluxm, 1, 0, myThid )
                0106       _EXCH_XY_RL(whfluxm   , myThid )
cf705a6c8e Mart*0107       DO bj = myByLo(myThid),myByHi(myThid)
                0108         DO bi = myBxLo(myThid),myBxHi(myThid)
                0109           DO j = 1-OLy,sNy+OLy
                0110             DO i = 1-OLx,sNx+OLx
0e09621e3e Patr*0111 c            print*,'Uncertainties for Heat Flux',i,j,whfluxm(i,j,bi,bj)
                0112              IF (whfluxm(i,j,bi,bj) .NE. 0. _d 0) THEN
                0113                  whfluxm(i,j,bi,bj) = 1. _d 0 /whfluxm(i,j,bi,bj)
                0114      &                                        /whfluxm(i,j,bi,bj)
                0115              ELSE
                0116                  whfluxm(i,j,bi,bj) = 1. _d 0
                0117              ENDIF
                0118             ENDDO
                0119           ENDDO
                0120         ENDDO
                0121       ENDDO
                0122 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                0123       CALL ACTIVE_WRITE_XY('whfluxm',whfluxm,1,0,myThid,dummy)
                0124 #endif
2e7aec9951 dngo*0125 #endif /* ALLOW_COST_HFLUXM */
0e09621e3e Patr*0126       RETURN
                0127       END