pkg/ecco/cost_gencost_transp.F

9c576a35a3 An T*0001 #include "ECCO_OPTIONS.h"
                0002
ebcc29af97 Jean*0003       SUBROUTINE COST_GENCOST_TRANSP( myThid )
9c576a35a3 An T*0004
ebcc29af97 Jean*0005 C     ==================================================================
                0006 C     SUBROUTINE COST_GENCOST_TRANSP
                0007 C     ==================================================================
                0008 C     o Evaluate cost function contributions of section transport.
9c576a35a3 An T*0009
ebcc29af97 Jean*0010       IMPLICIT NONE
9c576a35a3 An T*0011
ebcc29af97 Jean*0012 C     == global variables ==
9c576a35a3 An T*0013 #include "EEPARAMS.h"
                0014 #include "SIZE.h"
                0015 #include "PARAMS.h"
                0016 #include "GRID.h"
                0017 #ifdef ALLOW_CAL
                0018 # include "cal.h"
                0019 #endif
                0020 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0021 # include "ECCO_SIZE.h"
                0022 # include "ECCO.h"
9c576a35a3 An T*0023 #endif
                0024
ebcc29af97 Jean*0025 C     == routine arguments ==
                0026       INTEGER myThid
9c576a35a3 An T*0027
                0028 #ifdef ALLOW_GENCOST_CONTRIBUTION
ebcc29af97 Jean*0029 C     == external functions ==
                0030       INTEGER  ILNBLNK
                0031       EXTERNAL ILNBLNK
9c576a35a3 An T*0032
ebcc29af97 Jean*0033 C     == local variables ==
                0034       INTEGER nrecloc, localrec
                0035       INTEGER localstartdate(4)
9c576a35a3 An T*0036
ebcc29af97 Jean*0037       _RL myobs   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0038       _RL mybar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0039       _RL localdif(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0040       _RL difmask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0041       _RL localtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
9c576a35a3 An T*0042
ebcc29af97 Jean*0043       _RL dummy, facCost, facNum
9c576a35a3 An T*0044       _RL localperiod
                0045       _RL spminloc, spmaxloc, spzeroloc
ebcc29af97 Jean*0046       _RL tmpMeanTile(nSx,nSy), tmpNumTile(nSx,nSy)
                0047       _RL tmpMeanGlo, tmpNumGlo
                0048
                0049       INTEGER kgen(NGENCOST3D)
                0050       INTEGER bi, bj
                0051       INTEGER k
                0052       INTEGER obsrec, irec
                0053       INTEGER il,k2
                0054       INTEGER icount, icount_transp
                0055       LOGICAL dosumsq, dovarwei, doreadobs
                0056       LOGICAL exst
                0057
                0058       INTEGER preproc_i(NGENPPROC)
9c576a35a3 An T*0059       _RL preproc_r(NGENPPROC)
ebcc29af97 Jean*0060       CHARACTER*(MAX_LEN_FNAM) mybarfile
                0061       CHARACTER*(MAX_LEN_FNAM) preproc(NGENPPROC)
                0062       CHARACTER*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
                0063       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0064       CHARACTER*(128) fname1, fname0
                0065 C     == end of interface ==
9c576a35a3 An T*0066
ebcc29af97 Jean*0067 C=============== PART 0: initilization ===================
9c576a35a3 An T*0068
ebcc29af97 Jean*0069 C- facNum is 1 divided by the number of tiles in SIZE dot h
                0070       facNum = nSx*nPx
                0071       facNum = 1. _d 0 / facNum
9c576a35a3 An T*0072
ebcc29af97 Jean*0073 C-- detect the relevant gencost indices
8af4a14f0f An T*0074       do k=1,NGENCOST3D
                0075         kgen(k)=0
                0076       enddo
                0077
ebcc29af97 Jean*0078 C-- write a report of how many transport costs
8af4a14f0f An T*0079       write(msgbuf,'(A)') 'Inside cost_gencost_transp:'
                0080       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0081      &                    SQUEEZE_RIGHT, myThid )
                0082
                0083       icount_transp=0
9c576a35a3 An T*0084       do k=1,NGENCOST
                0085         if ( (gencost_name(k)(1:6).EQ.'transp').AND.
8af4a14f0f An T*0086      &     (using_gencost(k)) ) then
                0087           icount_transp=icount_transp+1
                0088           kgen(icount_transp)=k
ebcc29af97 Jean*0089           il=ILNBLNK(gencost_barfile(kgen(icount_transp)))
8af4a14f0f An T*0090           write(msgbuf,'(A,i4,A,A)') 'Cost ',kgen(icount_transp),
                0091      &    ': ',gencost_barfile(kgen(icount_transp))(1:il)
                0092           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0093      &                        SQUEEZE_RIGHT, myThid )
                0094
                0095         endif
9c576a35a3 An T*0096       enddo
                0097
ebcc29af97 Jean*0098 C-- write a report of how many transport costs
8af4a14f0f An T*0099
                0100       do icount=1,icount_transp
ebcc29af97 Jean*0101        if (kgen(icount).NE.0) then
9c576a35a3 An T*0102
ebcc29af97 Jean*0103 C ========
9c576a35a3 An T*0104
ebcc29af97 Jean*0105 C-- initialize objf and num:
                0106         DO bj = myByLo(myThid), myByHi(myThid)
                0107          DO bi = myBxLo(myThid), myBxHi(myThid)
8af4a14f0f An T*0108             objf_gencost(bi,bj,kgen(icount))=0. _d 0
                0109             num_gencost(bi,bj,kgen(icount))=0. _d 0
9c576a35a3 An T*0110          ENDDO
                0111         ENDDO
                0112
ebcc29af97 Jean*0113 C--   Initialise local variables.
9c576a35a3 An T*0114         nrecloc=0
8af4a14f0f An T*0115         nrecloc=gencost_nrec(kgen(icount))
9c576a35a3 An T*0116
ebcc29af97 Jean*0117 C-- only enters if there is at least 1 record
9c576a35a3 An T*0118         if(nrecloc.gt.0) then
                0119
ebcc29af97 Jean*0120           facCost = nrecloc
                0121           facCost = 1. _d 0 / facCost
9c576a35a3 An T*0122
                0123           localperiod=0.
8af4a14f0f An T*0124           localperiod=gencost_period(kgen(icount))
                0125           dummy=gencost_dummy(kgen(icount))
                0126           spminloc=gencost_spmin(kgen(icount))
                0127           spmaxloc=gencost_spmax(kgen(icount))
                0128           spzeroloc=gencost_spzero(kgen(icount))
                0129
ebcc29af97 Jean*0130 C prefer to have preproc match nosumsq but can not seem to get syntax
                0131 C to work for comparison of characters so match dosumsq for now.
9c576a35a3 An T*0132           dosumsq=.FALSE.
                0133           dovarwei=.FALSE.
                0134           do k2 = 1, NGENPPROC
8af4a14f0f An T*0135             preproc(k2)=gencost_preproc(k2,kgen(icount))
                0136             preproc_i(k2)=gencost_preproc_i(k2,kgen(icount))
                0137             preproc_c(k2)=gencost_preproc_c(k2,kgen(icount))
                0138             preproc_r(k2)=gencost_preproc_r(k2,kgen(icount))
9c576a35a3 An T*0139             if (preproc(k2).EQ.'variaweight') dovarwei=.TRUE.
8af4a14f0f An T*0140             if (preproc(k2)(1:7).EQ.'dosumsq') dosumsq=.TRUE.
9c576a35a3 An T*0141           enddo
                0142
ebcc29af97 Jean*0143 C-- report of dosumsq flag to make sure it is false
                0144           il=ILNBLNK(gencost_name(kgen(icount)))
                0145           write(msgBuf,'(3A,L5,2A)')
8af4a14f0f An T*0146      &    'Cost ',gencost_name(kgen(icount))(1:il),
                0147      &    ' dosumsq: ',dosumsq,' preproc(1): ',preproc(1)(1:7)
fce41e6d01 An T*0148           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0149      &                        SQUEEZE_RIGHT, myThid )
9c576a35a3 An T*0150
ebcc29af97 Jean*0151 C-- set mybarfile
8af4a14f0f An T*0152           mybarfile=gencost_barfile(kgen(icount))
                0153
ebcc29af97 Jean*0154 C-- set obsfile if defined. Not use for now. send warning
9c576a35a3 An T*0155           doreadobs=.FALSE.
8af4a14f0f An T*0156           if( .not. gencost_datafile(kgen(icount)).eq.' ') then
                0157 c            doreadobs=.TRUE.
ebcc29af97 Jean*0158             write(msgBuf,'(A)')
8af4a14f0f An T*0159      &      '**WARNING** S/R COST_GENCOST_TRANSP: gencost_datafile '
                0160             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0161      &                          SQUEEZE_RIGHT, myThid )
ebcc29af97 Jean*0162             write(msgBuf,'(A)')
8af4a14f0f An T*0163      &      'are currently ignored. Adjust the S/R to add code.'
                0164             CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0165      &                          SQUEEZE_RIGHT, myThid )
                0166           endif
9c576a35a3 An T*0167
ebcc29af97 Jean*0168 C model mask[W,S]: already included in transp calc in ecco_phys
9c576a35a3 An T*0169
ebcc29af97 Jean*0170 C--------- PART 0.1 read weights --------------------
8af4a14f0f An T*0171
ebcc29af97 Jean*0172 C-- read mask in gencost_mask
                0173 C-- for now, assume non-time-variable mask
9c576a35a3 An T*0174
ebcc29af97 Jean*0175 C=============== PART 1: main loop ===================
9c576a35a3 An T*0176           do irec = 1,nrecloc
8af4a14f0f An T*0177
ebcc29af97 Jean*0178 C--------- PART 1.1 read barfiles ------------------
8af4a14f0f An T*0179
ebcc29af97 Jean*0180 C-- set all bars to zeros:
fce41e6d01 An T*0181             call ecco_zero(mybar,Nr,zeroRL,myThid)
9c576a35a3 An T*0182
ebcc29af97 Jean*0183 C gencost_mask and fname0 are dummy, get fname1 from mybarfile
9c576a35a3 An T*0184             exst=.FALSE.
                0185             call ecco_zero(localtmp,Nr,zeroRL,myThid)
736e27304c Timo*0186             call cost_gencal(mybarfile,gencost_mask(kgen(icount)),
fce41e6d01 An T*0187      &       irec,localstartdate,localperiod,fname1,
ebcc29af97 Jean*0188      &       fname0,localrec,obsrec,exst,myThid)
                0189             call cost_genread(fname1,mybar,localtmp,irec,Nr,Nr,
9c576a35a3 An T*0190      &       nrecloc,preproc,preproc_c,preproc_i,preproc_r,
ebcc29af97 Jean*0191      &       dummy,myThid)
9c576a35a3 An T*0192
ebcc29af97 Jean*0193 C--------- PART 1.2 read data --------------------
8af4a14f0f An T*0194
ebcc29af97 Jean*0195 C-- ignore for now, but use doreadobs flag if needed
                0196 C-- be careful of recomputation when put inside if-end block
fce41e6d01 An T*0197 c            if(doreadobs) then
                0198             call ecco_zero(myobs,Nr,zeroRL,myThid)
                0199 c            endif
9c576a35a3 An T*0200
ebcc29af97 Jean*0201 C--------- PART 1.3 Cost calculation -------------
8af4a14f0f An T*0202
ebcc29af97 Jean*0203 C-- keep total at each irec to print out for time-series
                0204             DO bj = myByLo(myThid), myByHi(myThid)
                0205              DO bi = myBxLo(myThid), myBxHi(myThid)
                0206                tmpMeanTile(bi,bj) = 0. _d 0
                0207                tmpNumTile(bi,bj) = 0. _d 0
                0208              ENDDO
9c576a35a3 An T*0209             ENDDO
                0210
ebcc29af97 Jean*0211 C compute obs minus bar (localdif) and mask (difmask)
                0212 C note localtmp is set to 1.
                0213             call ecco_zero(localtmp,Nr,oneRL,myThid)
                0214             call ecco_zero(localdif,Nr,zeroRL,myThid)
                0215             call ecco_zero(difmask,Nr,zeroRL,myThid)
9c576a35a3 An T*0216
ebcc29af97 Jean*0217 C take care to set sp[min,max,zero]loc carefully to not
                0218 C filter out signal.  Can consider skip diffmsk step
9c576a35a3 An T*0219             call ecco_diffmsk(
11c3150c71 Mart*0220      I       mybar, myobs, localtmp,
                0221      I       Nr, Nr,spminloc, spmaxloc, spzeroloc,
fce41e6d01 An T*0222      O       localdif, difmask,
9c576a35a3 An T*0223      I       myThid )
                0224
ebcc29af97 Jean*0225             call ecco_zero(localtmp,Nr,oneRL,myThid)
9c576a35a3 An T*0226             call ecco_addcost(
ebcc29af97 Jean*0227      I       localdif,localtmp,difmask,Nr,Nr,dosumsq,
9c576a35a3 An T*0228      O       tmpMeanTile,tmpNumTile,
ebcc29af97 Jean*0229      I       myThid)
9c576a35a3 An T*0230
ebcc29af97 Jean*0231 cC either use ecco_diffmsk and ecco_addcost from above or simplify
                0232 cC to call below. For now keep syntax consistent with gencost
                0233 c            call ecco_zero(localtmp,Nr,oneRL,myThid)
8af4a14f0f An T*0234 c            call ecco_addcost(
ebcc29af97 Jean*0235 c     I       mybar,localtmp,localtmp,Nr,Nr,dosumsq,
8af4a14f0f An T*0236 c     O       tmpMeanTile,tmpNumTile,
ebcc29af97 Jean*0237 c     I       myThid)
8af4a14f0f An T*0238
ebcc29af97 Jean*0239 C global sums for display of time series
                0240 C note tmpNumGlo is the constant total wet points in the gencost_mask[W,S]
9c576a35a3 An T*0241             tmpMeanGlo = 0. _d 0
                0242             tmpNumGlo = 0. _d 0
ebcc29af97 Jean*0243             il=ILNBLNK(gencost_barfile(kgen(icount)))
9c576a35a3 An T*0244             CALL GLOBAL_SUM_TILE_RL( tmpMeanTile, tmpMeanGlo, myThid )
                0245             CALL GLOBAL_SUM_TILE_RL( tmpNumTile, tmpNumGlo, myThid )
8af4a14f0f An T*0246             WRITE(msgBuf,'(2A,I3,A,1PE21.14,1PE21.14)')
                0247      &        'globalsum transp ',gencost_barfile(kgen(icount))(1:il),
                0248      &        irec,' ',tmpMeanGlo,tmpNumGlo
9c576a35a3 An T*0249             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0250      &                          SQUEEZE_RIGHT, myThid )
                0251
ebcc29af97 Jean*0252 C sum that is actually be used in cost function
                0253             DO bj = myByLo(myThid), myByHi(myThid)
                0254              DO bi = myBxLo(myThid), myBxHi(myThid)
                0255                objf_gencost(bi,bj,kgen(icount))=
                0256      &            objf_gencost(bi,bj,kgen(icount))+tmpMeanTile(bi,bj)
                0257              ENDDO
9c576a35a3 An T*0258             ENDDO
                0259
                0260           enddo !irec
                0261
ebcc29af97 Jean*0262 C-- last step:
                0263 C-- divide by number of record to get mean transport:
                0264 C-- make num_gencost equals number of months/days used
                0265           DO bj = myByLo(myThid), myByHi(myThid)
                0266            DO bi = myBxLo(myThid), myBxHi(myThid)
                0267              objf_gencost(bi,bj,kgen(icount))=
                0268      &          objf_gencost(bi,bj,kgen(icount))*facCost
                0269              num_gencost(bi,bj,kgen(icount))=nrecloc*facNum
                0270            ENDDO
                0271           ENDDO
9c576a35a3 An T*0272
                0273         endif !if (nrecloc.gt.0)
ebcc29af97 Jean*0274        endif !if (kgen.NE.0)
8af4a14f0f An T*0275       enddo !icount_transp
9c576a35a3 An T*0276
                0277 #endif /* ALLOW_GENCOST_CONTRIBUTION */
                0278
                0279       RETURN
ebcc29af97 Jean*0280       END