Back to home page

MITgcm
 
 

    

File indexing completed on 2022-01-06 06:12:37 UTC

view on githubraw file Latest commit 9f5240b5 on 2022-01-05 15:24:45 UTC
df462307fb Timo*0001 #include "ECCO_OPTIONS.h"
                0002 
9f5240b52a Jean*0003       subroutine cost_gencost_moc(myThid)
df462307fb Timo*0004 
                0005 c     ==================================================================
                0006 c     SUBROUTINE cost_gencost_moc
                0007 c     ==================================================================
                0008 c
                0009 c     o Evaluate cost function contributions from MOC defined:
                0010 c
                0011 c       MOC_max = max_k { cumsum_k { zonally integrated meridional volume transport } }
                0012 c
                0013 c     o mybar loads trVol from barfile, which is masked by S/W edge mask
                0014 c       these masks denote the latitude line for taking zonal integral
                0015 c
                0016 c     ==================================================================
                0017 c     SUBROUTINE cost_gencost_moc
                0018 c     ==================================================================
                0019 
                0020       implicit none
                0021 
                0022 c     == global variables ==
                0023 
                0024 #include "EEPARAMS.h"
                0025 #include "SIZE.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
                0028 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0029 # include "ECCO_SIZE.h"
                0030 # include "ECCO.h"
df462307fb Timo*0031 #endif
                0032 
                0033 c     == routine arguments ==
9f5240b52a Jean*0034       integer myThid
df462307fb Timo*0035 
                0036 #ifdef ALLOW_GENCOST_CONTRIBUTION
                0037 
                0038 c     == local variables ==
                0039 
9f5240b52a Jean*0040       integer nrecloc, ioUnit
df462307fb Timo*0041 
9f5240b52a Jean*0042       _RL mybar     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
efbf0849bd Jean*0043       _RL gencost_mskTemporal
                0044       _RL tmpVar(1), dummyRL
df462307fb Timo*0045       _RS dummyRS(1)
9f5240b52a Jean*0046       _RL tmpCumSumTile(Nr,nSx,nSy)
df462307fb Timo*0047       _RL tmpNumTile(nSx,nSy)
9f5240b52a Jean*0048       _RL tmpCumSumGlo(Nr)
df462307fb Timo*0049       _RL tmpTile(nSx,nSy)
                0050       _RL myTempMax
                0051 
                0052       character*(MAX_LEN_FNAM) mybarfile
                0053       character*(MAX_LEN_MBUF) msgbuf
                0054       character*(128) fname0
                0055 
                0056       integer kgen, kg3
                0057       integer bi,bj
                0058       integer i,j,k
                0059       integer itlo,ithi
                0060       integer jtlo,jthi
                0061       integer imin, imax
                0062       integer jmin, jmax
                0063       integer irec
                0064       integer il
                0065       integer kmax
                0066       logical doglobalread
                0067       logical ladinit
                0068       logical exst
                0069 
                0070 c     == external functions ==
                0071 
                0072       integer  ilnblnk
                0073       external ilnblnk
                0074 
                0075 c     == end of interface ==
                0076 
9f5240b52a Jean*0077       jtlo = myByLo(myThid)
                0078       jthi = myByHi(myThid)
                0079       itlo = myBxLo(myThid)
                0080       ithi = myBxHi(myThid)
df462307fb Timo*0081       imin = 1
                0082       imax = sNx
                0083       jmin = 1
                0084       jmax = sNy
                0085 
                0086 CADJ  INIT tapelev_moc   = common, Nr
                0087 
                0088       do kgen=1,NGENCOST
                0089         kg3 = gencost_pointer3d(kgen)
                0090 
                0091 c-- detect the relevant gencost indices
                0092         if ( (gencost_name(kgen)(1:3).EQ.'moc').AND.
                0093      &     (using_gencost(kgen)) ) then
                0094 
                0095           write(msgbuf,'(A)') 'Inside cost_gencost_moc ...'
                0096           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0097      &                    SQUEEZE_RIGHT, myThid )
                0098 
                0099           il=ilnblnk(gencost_barfile(kgen))
                0100           write(msgbuf,'(A,i4,A,A)') 'Cost ',kgen,
                0101      &    ': ',gencost_barfile(kgen)(1:il)
                0102           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0103      &                        SQUEEZE_RIGHT, myThid )
                0104 
                0105 c=============== PART 0: initilization ===================
                0106 
                0107 c-- local variables
                0108 
                0109         doglobalread = .false.
                0110         ladinit = .false.
                0111         dummyRL = gencost_dummy(kgen)
                0112         nrecloc=gencost_nrec(kgen)
                0113         il=ilnblnk(gencost_barfile(kgen))
                0114         write(mybarfile,'(2a,i10.10)')
                0115      &    gencost_barfile(kgen)(1:il),'.',eccoiter
                0116 
                0117 c-- Initialize to zero
                0118           DO bj = jtlo,jthi
                0119             DO bi = itlo,ithi
                0120               do k = 1, Nr
                0121                 tmpCumSumTile(k,bi,bj) = 0. _d 0
                0122               enddo
                0123               tmpNumTile(bi,bj) = 0. _d 0
                0124             ENDDO
                0125           ENDDO
                0126 c model mask[W,S]: already included in transp calc in ecco_phys
                0127         nrecloc=gencost_nrec(kgen)
                0128 
                0129 c=============== PART 1: main loop ===================
                0130         do irec = 1,nrecloc
                0131 
                0132 c-- Read barfile
                0133 #ifdef ALLOW_AUTODIFF
                0134           call active_read_xyz( mybarfile, mybar, irec,
                0135      &                       doglobalread, ladinit,
9f5240b52a Jean*0136      &                       eccoiter, myThid,
df462307fb Timo*0137      &                       dummyRL )
                0138 #else
                0139           call READ_REC_XYZ_RL( mybarfile, mybar, irec,
9f5240b52a Jean*0140      &                       1, myThid )
df462307fb Timo*0141 #endif /* ALLOW_AUTODIFF */
                0142 
                0143 c-- Initialize after read
                0144 
                0145           DO bj = jtlo,jthi
                0146             DO bi = itlo,ithi
                0147               do k = 1, Nr
                0148                 tmpCumSumTile(k,bi,bj) = 0. _d 0
                0149               enddo
                0150               tmpNumTile(bi,bj) = 0. _d 0
                0151             ENDDO
                0152           ENDDO
                0153 
                0154 c-- Temporal mask
                0155           il = ilnblnk(gencost_mask(kgen))
                0156           write(fname0(1:128),'(2A)')
                0157      &      gencost_mask(kgen)(1:il),'T'
                0158           inquire( file=fname0(1:il+1), exist=exst )
                0159 
                0160           if ( (.NOT.exst).OR.(gencost_mask(kgen).EQ.' ')
                0161      &          ) then
                0162 
                0163            write(msgBuf,'(3A)') '**Warning: temporal msk file: ',
                0164      &       fname0(1:il+1), ' not found, using 1/nrecloc'
                0165            CALL PRINT_MESSAGE(msgBuf, standardMessageUnit,
                0166      &         SQUEEZE_RIGHT, myThid )
                0167 
                0168            gencost_mskTemporal=nrecloc
                0169            gencost_mskTemporal=1. _d 0 / gencost_mskTemporal
                0170           else
                0171 
                0172            write(msgBuf,'(2A)') 'Using temporal msk from file: ',
                0173      &          fname0(1:il+1)
                0174            CALL PRINT_MESSAGE(msgBuf, standardMessageUnit,
                0175      &          SQUEEZE_RIGHT, myThid )
                0176 
                0177            ioUnit = 0
efbf0849bd Jean*0178            CALL MDS_READVEC_LOC( fname0, cost_iprec, ioUnit, 'RL',
                0179      &                      1, tmpVar, dummyRS, 0, 0, irec, myThid )
                0180            gencost_mskTemporal = tmpVar(1)
df462307fb Timo*0181           endif
                0182 
                0183 c=============== PART 2: Cost Computation ===================
                0184 c-- Compute cost only if nonzero temporal mask
                0185           if ( gencost_mskTemporal .ne. 0 ) then
                0186 
                0187           if ( myProcId .EQ. 0 ) num_gencost(1,1,kgen)=
                0188      &      num_gencost(1,1,kgen)+gencost_mskTemporal
                0189 
                0190 c=============== PART 2.1: Cumulative sum ===================
                0191 c Take cumulative sum of my bar from bottom up
                0192 c i.e. compute the streamfunction (assuming mybar = trvol)
                0193 c
                0194           DO bj = jtlo,jthi
                0195             DO bi = itlo,ithi
                0196               do k = Nr, 1, -1
                0197                 do j = jmin,jmax
                0198                   do i = imin,imax
                0199                     tmpCumSumTile(k,bi,bj)=tmpCumSumTile(k,bi,bj) -
                0200      &                mybar(i,j,k,bi,bj)*gencost_mskTemporal
                0201 
                0202                     tmpNumTile(bi,bj)=
                0203      &                tmpNumTile(bi,bj)+1. _d 0
                0204                   enddo
                0205                 enddo
                0206                if( k .lt. Nr ) then
                0207                 tmpCumSumTile(k,bi,bj)=tmpCumSumTile(k,bi,bj) +
                0208      &            tmpCumSumTile(k+1,bi,bj)
                0209                endif
                0210               enddo
                0211             ENDDO
                0212           ENDDO
                0213 
                0214 c-- Compute global sum at each level
                0215           do k = Nr, 1, -1
                0216             DO bj = jtlo,jthi
                0217               DO bi = itlo,ithi
                0218                 tmpTile(bi,bj) = tmpCumSumTile(k,bi,bj)
                0219               ENDDO
                0220             ENDDO
                0221             tmpCumSumGlo(k) = 0. _d 0
9f5240b52a Jean*0222             CALL GLOBAL_SUM_TILE_RL(tmpTile, tmpCumSumGlo(k), myThid )
df462307fb Timo*0223           enddo
                0224 
                0225 c=============== PART 2.2: Get max val ===================
                0226 c-- Find maximum in global cumulative sum
                0227 
                0228           myTempMax = tmpCumSumGlo(1)
                0229           kmax = 1
                0230 
                0231           do k = 2, Nr
                0232 CADJ STORE myTempMax = tapelev_moc, key = k
                0233             if( myTempMax < tmpCumSumGlo(k) ) then
                0234               myTempMax = tmpCumSumGlo(k)
                0235               kmax = k
                0236             endif
                0237           enddo
                0238 
                0239           WRITE(msgBuf,'(2A,I3,A,1PE21.14,A,I2)')
                0240      &        'moc cost ',gencost_barfile(kgen)(1:il),
                0241      &        irec,' ', myTempMax, 'kmax: ',kmax
                0242           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0243      &                          SQUEEZE_RIGHT, myThid )
                0244 
                0245 c=============== PART 2.3: Save to obj ===================
                0246 c-- Add MOC contribution to actual objective function variable
                0247 c   Note: using global values, so only want to store in
                0248 c         one processors obj function value
                0249 
                0250           if ( myProcId .EQ. 0 ) objf_gencost(1,1,kgen)=
                0251      &      objf_gencost(1,1,kgen)+myTempMax
                0252 
                0253 c============= Done with cost computation =====================
                0254 
                0255           else ! mskTemporal == 0
                0256 
                0257             WRITE(msgBuf,'(A,I3,A,I3)')
                0258      &          'gencost_mskTemporal = 0, irec: ',irec, ' / ',
                0259      &          nrecloc
                0260           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0261      &                          SQUEEZE_RIGHT, myThid )
                0262 
                0263           endif ! mskTemporal /=0
                0264           enddo ! irec=1->nrecloc
                0265 
                0266 c-- Print out what actually is used as cost function
                0267           WRITE(msgBuf,'(A,1PE21.14)') 'moc fc: ',
                0268      &          objf_gencost(1,1,kgen)
                0269           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0270      &                        SQUEEZE_RIGHT, myThid )
                0271 
                0272         endif ! gencost_name(kgen)=moc
                0273       enddo ! kgen=1->NGENCOST
                0274 
                0275 #endif /* ALLOW_GENCOST_CONTRIBUTION */
                0276 
                0277       return
                0278       end
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!