Back to home page

MITgcm
 
 

    

File indexing completed on 2023-07-14 05:10:22 UTC

view on githubraw file Latest commit de57a2ec on 2023-07-13 16:55:13 UTC
f09238ab8f Gael*0001 #include "ECCO_OPTIONS.h"
                0002 
                0003       subroutine cost_gencost_bpv4(
4890e8ebab antn*0004      I                     myThid )
f09238ab8f Gael*0005 
                0006 c     ==================================================================
                0007 c     SUBROUTINE cost_gencost_bpv4
                0008 c     ==================================================================
                0009 c
                0010 c     o Evaluate cost function contribution of bottom pressure anoamlies
                0011 c       => GRACE data
                0012 c
                0013 c     started: Gael Forget Oct-2009
                0014 c
                0015 c     ==================================================================
                0016 c     SUBROUTINE cost_bp
                0017 c     ==================================================================
                0018 
                0019       implicit none
                0020 
                0021 c     == global variables ==
                0022 
                0023 #include "EEPARAMS.h"
                0024 #include "SIZE.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "DYNVARS.h"
                0028 
                0029 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0030 # include "ECCO_SIZE.h"
                0031 # include "ECCO.h"
f09238ab8f Gael*0032 #endif
                0033 
                0034 c     == routine arguments ==
f8e779c983 antn*0035       integer myThid
f09238ab8f Gael*0036 
                0037 #ifdef ALLOW_ECCO
                0038 #ifdef ALLOW_GENCOST_CONTRIBUTION
                0039 
4890e8ebab antn*0040 c     == external functions ==
                0041       integer  ilnblnk
                0042       external ilnblnk
f09238ab8f Gael*0043 
4890e8ebab antn*0044 c     == local variables ==
f09238ab8f Gael*0045       integer bi,bj
                0046       integer i,j
                0047       integer irec
                0048       integer il
                0049       logical doglobalread
                0050       logical ladinit
f8e779c983 antn*0051       _RL locbpbar(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0052       _RL locbpdat(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0053       _RL locbpmask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0054       _RL locwbp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0055       _RL bpdifmean ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
                0056       _RL bpdifanom ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
                0057       _RL bpdatmean ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
                0058       _RL bpdatanom ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
                0059       _RL bpcount ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
f09238ab8f Gael*0060       _RL junk
de57a2ec4b Mart*0061       character*(MAX_LEN_FNAM) fname
                0062       character*(MAX_LEN_FNAM) fname4test
f09238ab8f Gael*0063       _RL fac
                0064       _RL offset
                0065       _RL offset_sum
                0066       integer k, kgen
b9b1c4d04c Timo*0067       logical dosumsq
4890e8ebab antn*0068 #ifdef ALLOW_SMOOTH
                0069       integer k2_bp
                0070       logical exst
                0071 #endif
                0072 #ifdef ECCO_VERBOSE
                0073       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0074 #endif
f09238ab8f Gael*0075 c     == end of interface ==
                0076 
                0077       kgen=0
                0078       do k=1,NGENCOST
b304983705 Gael*0079         if ( (gencost_name(k).EQ.'bpv4-grace').AND.
6b2230d510 Ou W*0080      &       (.NOT.gencost_is1d(k)).AND.
b304983705 Gael*0081      &       (using_gencost(k)) ) kgen=k
f09238ab8f Gael*0082       enddo
                0083 
                0084       if (kgen.GT.0) then
                0085 
4890e8ebab antn*0086 #ifdef ALLOW_SMOOTH
                0087        k2_bp=0
                0088        if ( useSMOOTH) then
                0089           do k = 1,ngenpproc
                0090             if (gencost_posproc(k,kgen).eq.'smooth') k2_bp=k
                0091           enddo
                0092        endif
                0093 #endif
                0094 
b9b1c4d04c Timo*0095        dosumsq=.TRUE.
11c3150c71 Mart*0096        call ecco_zero( gencost_weight(1-OLx,1-OLy,1,1,kgen),
                0097      &                 1, zeroRL, myThid )
f8e779c983 antn*0098        if ( gencost_errfile(kgen) .NE. ' ' )
11c3150c71 Mart*0099      &   call ecco_readwei( gencost_errfile(kgen),
                0100      &     gencost_weight(1-OLx,1-OLy,1,1,kgen),
                0101      &     1, 1, 1, dosumsq, myThid )
cbd85e4123 Gael*0102 
f09238ab8f Gael*0103 c-- initialise local variables
                0104 cgf convert phibot from m2/s2 to cm
f8e779c983 antn*0105        fac = 1. _d 2 / 9.81 _d 0
4890e8ebab antn*0106        DO bj = myByLo(myThid), myByHi(myThid)
                0107         DO bi = myBxLo(myThid), myBxHi(myThid)
f8e779c983 antn*0108           do j = 1,sNy
                0109            do i = 1,sNx
                0110              bpdifmean(i,j,bi,bj) = 0. _d 0
                0111              bpdifanom(i,j,bi,bj) = 0. _d 0
                0112              bpdatmean(i,j,bi,bj) = 0. _d 0
                0113              bpdatanom(i,j,bi,bj) = 0. _d 0
                0114              bpcount(i,j,bi,bj) = 0. _d 0
                0115              locwbp(i,j,bi,bj) = 0. _d 0
                0116              locbpbar(i,j,bi,bj) = 0. _d 0
                0117              locbpdat(i,j,bi,bj) = 0. _d 0
                0118              locbpmask(i,j,bi,bj) = 0. _d 0
                0119            enddo
f09238ab8f Gael*0120           enddo
4890e8ebab antn*0121         ENDDO
                0122        ENDDO
f09238ab8f Gael*0123 
f8e779c983 antn*0124        doglobalread = .false.
                0125        ladinit      = .false.
f09238ab8f Gael*0126 
                0127 c-- map global variable to local variables
                0128 
4890e8ebab antn*0129        DO bj = myByLo(myThid), myByHi(myThid)
                0130         DO bi = myBxLo(myThid), myBxHi(myThid)
f8e779c983 antn*0131           do j = 1,sNy
                0132            do i = 1,sNx
                0133              locwbp(i,j,bi,bj) = gencost_weight(i,j,bi,bj,kgen)
                0134            enddo
                0135           enddo
4890e8ebab antn*0136         ENDDO
                0137        ENDDO
f09238ab8f Gael*0138 
                0139 #ifdef ALLOW_CTRL
f8e779c983 antn*0140        il=ilnblnk( gencost_barfile(kgen) )
de57a2ec4b Mart*0141        write(fname,'(2a,i10.10)')
f09238ab8f Gael*0142      &     gencost_barfile(kgen)(1:il),'.',eccoiter
                0143 #endif
                0144 
                0145 c--   ============
                0146 c--   Mean values.
                0147 c--   ============
                0148 
4890e8ebab antn*0149 #ifdef ECCO_VERBOSE
                0150        WRITE(msgBuf,'(A,1x,i5,1x,i10)')
                0151      &      ' bpv4, kgen, nmonsrec: ',kgen, nmonsrec
                0152        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0153      &                     SQUEEZE_RIGHT, myThid )
                0154 #endif
f8e779c983 antn*0155        do irec = 1, nmonsrec
f09238ab8f Gael*0156 
                0157 c--     Compute the mean over all bpdat records.
101f75e5cd Gael*0158 #ifdef ALLOW_AUTODIFF
f09238ab8f Gael*0159         call active_read_xy( fname, locbpbar, irec, doglobalread,
f8e779c983 antn*0160      &                       ladinit, eccoiter, myThid,
f09238ab8f Gael*0161      &                       gencost_dummy(kgen) )
101f75e5cd Gael*0162 #else
                0163         CALL READ_REC_XY_RL( fname, locbpbar,
                0164      &                       iRec, 1, myThid )
                0165 #endif
                0166 
f8e779c983 antn*0167         call cost_bp_read( gencost_datafile(kgen),
f09238ab8f Gael*0168      &       gencost_startdate(1,kgen),
f8e779c983 antn*0169      &       locbpdat, locbpmask, irec, myThid )
f09238ab8f Gael*0170 
4890e8ebab antn*0171 #ifdef ECCO_VERBOSE
                0172         il=ilnblnk( gencost_datafile(kgen) )
                0173         WRITE(msgBuf,'(A,1x,A,1x,i10)') ' bpv4, datafile, irec: ',
                0174      &        gencost_datafile(kgen)(1:il), irec
                0175         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0176      &                      SQUEEZE_RIGHT, myThid )
de57a2ec4b Mart*0177         write(fname4test,'(2a,i4.4)')
4890e8ebab antn*0178      &       gencost_datafile(kgen)(1:il),'.',irec
                0179         call write_rec_xy_rl( fname4test, locbpdat, 1, 1, myThid)
                0180 #endif
                0181 
                0182         DO bj = myByLo(myThid), myByHi(myThid)
                0183          DO bi = myBxLo(myThid), myBxHi(myThid)
f8e779c983 antn*0184            do j = 1,sNy
                0185             do i = 1,sNx
                0186               if ( (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
                0187      &             (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
                0188                 bpdifmean(i,j,bi,bj) = bpdifmean(i,j,bi,bj) +
f09238ab8f Gael*0189      &              ( fac*locbpbar(i,j,bi,bj) - locbpdat(i,j,bi,bj) )
f8e779c983 antn*0190                 bpdatmean(i,j,bi,bj) = bpdatmean(i,j,bi,bj) +
f09238ab8f Gael*0191      &              locbpdat(i,j,bi,bj)
f8e779c983 antn*0192                 bpcount(i,j,bi,bj) = bpcount(i,j,bi,bj) + 1. _d 0
                0193               endif
f09238ab8f Gael*0194             enddo
f8e779c983 antn*0195            enddo
4890e8ebab antn*0196          ENDDO
                0197         ENDDO
f09238ab8f Gael*0198 
f8e779c983 antn*0199        enddo
f09238ab8f Gael*0200 
4890e8ebab antn*0201        DO bj = myByLo(myThid), myByHi(myThid)
                0202         DO bi = myBxLo(myThid), myBxHi(myThid)
f8e779c983 antn*0203           do j = 1,sNy
                0204            do i = 1,sNx
                0205              if (bpcount(i,j,bi,bj).GT. 0. _d 0) then
                0206                bpdifmean(i,j,bi,bj) =
f09238ab8f Gael*0207      &              bpdifmean(i,j,bi,bj)/bpcount(i,j,bi,bj)
f8e779c983 antn*0208                bpdatmean(i,j,bi,bj) =
f09238ab8f Gael*0209      &              bpdatmean(i,j,bi,bj)/bpcount(i,j,bi,bj)
f8e779c983 antn*0210              endif
                0211            enddo
f09238ab8f Gael*0212           enddo
4890e8ebab antn*0213         ENDDO
                0214        ENDDO
f09238ab8f Gael*0215 
                0216 c--   ==========
                0217 c--   Anomalies.
                0218 c--   ==========
                0219 
                0220 c--   Loop over records for the second time.
f8e779c983 antn*0221        do irec = 1, nmonsrec
101f75e5cd Gael*0222 #ifdef ALLOW_AUTODIFF
f09238ab8f Gael*0223         call active_read_xy( fname, locbpbar, irec, doglobalread,
f8e779c983 antn*0224      &                       ladinit, eccoiter, myThid,
f09238ab8f Gael*0225      &                       gencost_dummy(kgen) )
101f75e5cd Gael*0226 #else
                0227         CALL READ_REC_XY_RL( fname, locbpbar,
                0228      &                       iRec, 1, myThid )
                0229 #endif
f09238ab8f Gael*0230 
f8e779c983 antn*0231         call cost_bp_read( gencost_datafile(kgen),
f09238ab8f Gael*0232      &       gencost_startdate(1,kgen),
f8e779c983 antn*0233      &       locbpdat, locbpmask, irec, myThid )
f09238ab8f Gael*0234 
                0235 c--    Compute field of anomalies
4890e8ebab antn*0236         DO bj = myByLo(myThid), myByHi(myThid)
                0237          DO bi = myBxLo(myThid), myBxHi(myThid)
f8e779c983 antn*0238            do j = 1,sNy
                0239             do i = 1,sNx
                0240               if ( (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
                0241      &             (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
                0242                 bpdifanom(i,j,bi,bj) =
f09238ab8f Gael*0243      &              ( fac*locbpbar(i,j,bi,bj) - locbpdat(i,j,bi,bj) )
                0244      &              - bpdifmean(i,j,bi,bj)
f8e779c983 antn*0245                 bpdatanom(i,j,bi,bj) =
f09238ab8f Gael*0246      &              locbpdat(i,j,bi,bj) - bpdatmean(i,j,bi,bj)
f8e779c983 antn*0247               else
                0248                 bpdifanom(i,j,bi,bj) = 0. _d 0
                0249                 bpdatanom(i,j,bi,bj) = 0. _d 0
                0250               endif
f09238ab8f Gael*0251             enddo
f8e779c983 antn*0252            enddo
4890e8ebab antn*0253          ENDDO
                0254         ENDDO
f09238ab8f Gael*0255 
                0256 c--    Remove global mean value
f8e779c983 antn*0257         offset     = 0. _d 0
                0258         offset_sum = 0. _d 0
f09238ab8f Gael*0259 
4890e8ebab antn*0260         DO bj = myByLo(myThid), myByHi(myThid)
                0261          DO bi = myBxLo(myThid), myBxHi(myThid)
f8e779c983 antn*0262            do j = 1,sNy
                0263             do i = 1,sNx
f09238ab8f Gael*0264               if ( (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
                0265      &             (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
                0266                 offset  = offset + RA(i,j,bi,bj)*bpdifanom(i,j,bi,bj)
                0267                 offset_sum = offset_sum + RA(i,j,bi,bj)
                0268               endif
                0269             enddo
f8e779c983 antn*0270            enddo
4890e8ebab antn*0271          ENDDO
                0272         ENDDO
f09238ab8f Gael*0273 
f8e779c983 antn*0274         _GLOBAL_SUM_RL( offset     , myThid )
                0275         _GLOBAL_SUM_RL( offset_sum , myThid )
f09238ab8f Gael*0276 
4890e8ebab antn*0277         DO bj = myByLo(myThid), myByHi(myThid)
                0278          DO bi = myBxLo(myThid), myBxHi(myThid)
f8e779c983 antn*0279            do j = 1,sNy
                0280             do i = 1,sNx
                0281               if ( (offset_sum.GT. 0. _d 0).AND.
                0282      &             (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
                0283      &             (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
                0284                 bpdifanom(i,j,bi,bj) = bpdifanom(i,j,bi,bj)
                0285      &                               - offset/offset_sum
                0286               endif
f09238ab8f Gael*0287             enddo
f8e779c983 antn*0288            enddo
4890e8ebab antn*0289          ENDDO
                0290         ENDDO
f09238ab8f Gael*0291 
                0292 c--    Smooth field of anomalies
f8e779c983 antn*0293         if (gencost_outputlevel(kgen).GT.0) then
de57a2ec4b Mart*0294          write(fname4test,'(1a)') 'bpdifanom_raw'
f8e779c983 antn*0295          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0296      &                         bpdifanom, irec, 1, myThid )
de57a2ec4b Mart*0297          write(fname4test,'(1a)') 'bpdatanom_raw'
f8e779c983 antn*0298          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0299      &                         bpdatanom, irec, 1, myThid )
                0300         endif
f09238ab8f Gael*0301 
                0302 #ifdef ALLOW_SMOOTH
4890e8ebab antn*0303         IF ( k2_bp.GE.1 ) THEN
                0304 C- Note: Unlike all other generic-cost smoothing, original code was using
                0305 C  "smooth_basic2D" with hard-coded smoothing length (3.e+5) and time
                0306 C  multiplicator (3000). Keep original call (now using run-time params setting)
                0307 C  when the scaling file gencost_posproc_c is missing:
                0308          il=ilnblnk( gencost_posproc_c(k2_bp,kgen) )
                0309          inquire(file=gencost_posproc_c(k2_bp,kgen)(1:il), exist=exst)
                0310          IF ( il .gt. 1 .and. exst ) THEN
                0311            call smooth_hetero2d(bpdifanom,maskc,
                0312      &      gencost_posproc_c(k2_bp,kgen),
                0313      &      gencost_posproc_i(k2_bp,kgen),myThid)
                0314          ELSE
                0315            call smooth_basic2D(bpdifanom,maskInC,
                0316      &      gencost_posproc_r(k2_bp,kgen),
                0317      &      gencost_posproc_i(k2_bp,kgen),myThid)
                0318          ENDIF
                0319         ENDIF
f09238ab8f Gael*0320 #endif
                0321 
f8e779c983 antn*0322         if (gencost_outputlevel(kgen).GT.0) then
f09238ab8f Gael*0323 #ifdef ALLOW_SMOOTH
4890e8ebab antn*0324          IF ( k2_bp.GE.1 ) THEN
                0325 C--   Skip second evaluation of "il" & "exst" (unchanged since first one above)
                0326 c         il=ilnblnk( gencost_posproc_c(k2_bp,kgen) )
                0327 c         inquire(file=gencost_posproc_c(k2_bp,kgen)(1:il), exist=exst)
                0328           IF ( il .gt. 1 .and. exst ) THEN
                0329            call smooth_hetero2d(bpdatanom,maskc,
                0330      &      gencost_posproc_c(k2_bp,kgen),
                0331      &      gencost_posproc_i(k2_bp,kgen),myThid)
                0332           ELSE
                0333            call smooth_basic2D(bpdatanom,maskInC,
                0334      &      gencost_posproc_r(k2_bp,kgen),
                0335      &      gencost_posproc_i(k2_bp,kgen),myThid)
                0336           ENDIF
                0337          ENDIF
f09238ab8f Gael*0338 #endif
                0339 
de57a2ec4b Mart*0340          write(fname4test,'(1a)') 'bpdifanom_smooth'
f8e779c983 antn*0341          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0342      &                         bpdifanom, irec, 1, myThid )
de57a2ec4b Mart*0343          write(fname4test,'(1a)') 'bpdatanom_smooth'
f8e779c983 antn*0344          CALL WRITE_REC_3D_RL( fname4test, precFloat32, 1,
                0345      &                         bpdatanom, irec, 1, myThid )
                0346         endif
f09238ab8f Gael*0347 
                0348 c--    Compute cost function
4890e8ebab antn*0349         DO bj = myByLo(myThid), myByHi(myThid)
                0350          DO bi = myBxLo(myThid), myBxHi(myThid)
f8e779c983 antn*0351            do j = 1,sNy
                0352             do i = 1,sNx
f09238ab8f Gael*0353 c-- map to global cost variables
f8e779c983 antn*0354               if ( (locwbp(i,j,bi,bj).NE. 0. _d 0).AND.
                0355      &             (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
                0356      &             (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
                0357                  junk = bpdifanom(i,j,bi,bj)
f4cfad5b95 Ou W*0358                  objf_gencost(bi,bj,kgen) = objf_gencost(bi,bj,kgen)
f09238ab8f Gael*0359      &               + junk*junk*locwbp(i,j,bi,bj)
f4cfad5b95 Ou W*0360                  num_gencost(bi,bj,kgen) = num_gencost(bi,bj,kgen)
f09238ab8f Gael*0361      &               + 1. _d 0
f8e779c983 antn*0362               endif
f09238ab8f Gael*0363             enddo
f8e779c983 antn*0364            enddo
4890e8ebab antn*0365          ENDDO
                0366         ENDDO
f09238ab8f Gael*0367 
f8e779c983 antn*0368        enddo
f09238ab8f Gael*0369 
                0370       endif !if (kgen.GT.0) then
                0371 
                0372 #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
                0373 #endif /* ifdef ALLOW_ECCO */
                0374 
f8e779c983 antn*0375       RETURN
                0376       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!