pkg/smooth/smooth_filtervar3d.F

0c3d35c9cd Gael*0001 #include "SMOOTH_OPTIONS.h"
                0002
9f5240b52a Jean*0003       SUBROUTINE SMOOTH_FILTERVAR3D( smoothOpNb, myThid )
0c3d35c9cd Gael*0004
                0005 C     *==========================================================*
                0006 C     | SUBROUTINE smooth_filtervar3D
                0007 C     | o Routine that computes the filter variance
                0008 C     |   field associated with a diffusion operator, as part
                0009 C     |   a 3D spatial correlation operator (smooth_correld3D.F)
                0010 C     |   See Weaver and Courtier 01 for details.
                0011 C     *==========================================================*
                0012
                0013       IMPLICIT NONE
9f5240b52a Jean*0014 c     == global variables ==
0c3d35c9cd Gael*0015 #include "SIZE.h"
                0016 #include "EEPARAMS.h"
                0017 #include "PARAMS.h"
                0018 #include "GRID.h"
                0019 #include "SMOOTH.h"
                0020
7b8b86ab99 Timo*0021       INTEGER smoothOpNb, myThid
5fda927278 Gael*0022
7b8b86ab99 Timo*0023 c     == external functions ==
                0024       REAL*8   port_rand, port_rand_norm
0410c4d6ff Jean*0025       EXTERNAL PORT_RAND, PORT_RAND_NORM
                0026
7b8b86ab99 Timo*0027 c     == local variables ==
                0028       INTEGER i,j,k, bi, bj, ii, jj, kk
                0029       INTEGER diLoc, djLoc,  dkLoc
                0030       INTEGER nbRand, nbt_in
de57a2ec4b Mart*0031       CHARACTER*(MAX_LEN_FNAM) fnamegeneric
0c3d35c9cd Gael*0032       _RL smoothTmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0033       _RL smoothTmpMean(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
7b8b86ab99 Timo*0034       _RL smoothTmpVar (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0035       _RS smooth3Dmask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
f9d7cbfb72 Ou W*0036
                0037       INTEGER IL
                0038
                0039 c     == functions ==
                0040       INTEGER  ILNBLNK
                0041       EXTERNAL ILNBLNK
                0042
7b8b86ab99 Timo*0043 c     == end of interface ==
0c3d35c9cd Gael*0044
f9d7cbfb72 Ou W*0045       IL = ILNBLNK( smoothDir )
                0046
7b8b86ab99 Timo*0047 c--   allow a mask other than maskC
9f5240b52a Jean*0048       DO bj=myByLo(myThid),myByHi(myThid)
                0049        DO bi=myBxLo(myThid),myBxHi(myThid)
7b8b86ab99 Timo*0050         DO k = 1,Nr
                0051          DO j = 1-OLy,sNy+OLy
                0052           DO i = 1-OLx,sNx+OLx
                0053            IF (smooth3DmaskName(smoothOpNb)(1:5).EQ.'maskC') THEN
                0054             smooth3Dmask(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
                0055            ELSEIF (smooth3DmaskName(smoothOpNb)(1:5).EQ.'maskW') THEN
                0056             smooth3Dmask(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)
                0057            ELSEIF (smooth3DmaskName(smoothOpNb)(1:5).EQ.'maskS') THEN
                0058             smooth3Dmask(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)
                0059            ENDIF
                0060           ENDDO
                0061          ENDDO
                0062         ENDDO
                0063        ENDDO
                0064       ENDDO
                0065
0c3d35c9cd Gael*0066 c if smooth3Dfilter(smoothOpNb)=0: the filter variance field
                0067 c has been computed earlier and is already in the run directory
                0068 c so this routine does not do anything
                0069
7b8b86ab99 Timo*0070       IF (smooth3Dfilter(smoothOpNb).NE.0) THEN
0c3d35c9cd Gael*0071
9f5240b52a Jean*0072        nbt_in = smooth3Dnbt(smoothOpNb)/2
0c3d35c9cd Gael*0073
                0074 c read smoothing [i.e diffusion] operator:
f9d7cbfb72 Ou W*0075        WRITE(fnamegeneric,'(2A,I3.3)')
                0076      &       smoothDir(1:IL),
9f5240b52a Jean*0077      &       'smooth3Doperator', smoothOpNb
                0078        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kwx,
                0079      &                      1, 1, myThid )
                0080        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kwy,
                0081      &                      2, 1, myThid )
                0082        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kwz,
                0083      &                      3, 1, myThid )
                0084        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kux,
                0085      &                      4, 1, myThid )
                0086        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kvy,
                0087      &                      5, 1, myThid )
                0088        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kuz,
                0089      &                      6, 1, myThid )
                0090        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kvz,
                0091      &                      7, 1, myThid )
                0092        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kuy,
                0093      &                      8, 1, myThid )
                0094        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr, smooth3D_Kvx,
                0095      &                      9, 1, myThid )
                0096        CALL READ_REC_3D_RL( fnamegeneric, smoothprec, Nr,
                0097      &                      smooth3D_kappaR, 10, 1, myThid )
                0098        CALL EXCH_XYZ_RL( smooth3D_Kwx, myThid )
                0099        CALL EXCH_XYZ_RL( smooth3D_Kwy, myThid )
                0100        CALL EXCH_XYZ_RL( smooth3D_Kwz, myThid )
                0101        CALL EXCH_XYZ_RL( smooth3D_Kux, myThid )
                0102        CALL EXCH_XYZ_RL( smooth3D_Kvy, myThid )
                0103        CALL EXCH_XYZ_RL( smooth3D_Kuz, myThid )
                0104        CALL EXCH_XYZ_RL( smooth3D_Kvz, myThid )
                0105        CALL EXCH_XYZ_RL( smooth3D_Kuy, myThid )
                0106        CALL EXCH_XYZ_RL( smooth3D_Kvx, myThid )
                0107        CALL EXCH_XYZ_RL( smooth3D_kappaR, myThid )
0c3d35c9cd Gael*0108
                0109 c initialize filter variance field:
9f5240b52a Jean*0110        DO bj=myByLo(myThid),myByHi(myThid)
                0111         DO bi=myBxLo(myThid),myBxHi(myThid)
                0112          DO k=1,Nr
                0113           DO j=1-OLy,sNy+OLy
                0114            DO i=1-OLx,sNx+OLx
                0115             smooth3Dnorm(i,j,k,bi,bj) = 0.
                0116            ENDDO
0c3d35c9cd Gael*0117           ENDDO
                0118          ENDDO
                0119         ENDDO
                0120        ENDDO
                0121
9f5240b52a Jean*0122        IF (smooth3Dfilter(smoothOpNb).EQ.2) THEN
0c3d35c9cd Gael*0123 c compute the normalization matrix using the approximate method
                0124 c
                0125 c This method can be quite expensive -- so that the approximate
                0126 c method (see below) is usually the prefered one.
                0127 c The exact method can be used to check the accuracy
                0128 c of the approximate method results (that can be predicted).
                0129 c
                0130 c note: the exact method requires the adjoint of smooth_diff2D.F (see below)
                0131
9f5240b52a Jean*0132         diLoc = 15 !int(5*smooth_L/smooth_dx)
                0133         djLoc = 20 !int(5*smooth_L/smooth_dx)
                0134         dkLoc = 8
                0135
                0136         DO kk=1,dkLoc
                0137          DO ii=1,diLoc,2
                0138           DO jj=1,djLoc,2
                0139
                0140            DO bj=myByLo(myThid),myByHi(myThid)
                0141             DO bi=myBxLo(myThid),myBxHi(myThid)
                0142              DO k=1,Nr
                0143               DO j=1-OLy,sNy+OLy
                0144                DO i=1-OLx,sNx+OLx
                0145                 smoothTmpFld(i,j,k,bi,bj) = 0.
                0146                ENDDO
                0147               ENDDO
                0148              ENDDO
                0149              DO k=kk,Nr,dkLoc
                0150               DO j=jj,sNy,djLoc
                0151                DO i=ii,sNx,diLoc
                0152                 smoothTmpFld(i,j,k,bi,bj) = 1.
                0153                ENDDO
                0154               ENDDO
                0155              ENDDO
                0156             ENDDO
                0157            ENDDO
0c3d35c9cd Gael*0158
                0159 c note: as we go to adjoint part, we need to have 0 in overlaps
                0160 c       so we must NOT have done an exchange for smoothTmpFld
                0161
                0162 c adjoint:
9f5240b52a Jean*0163            WRITE(errorMessageUnit,'(A,/,A)' )
                0164      &      'you need to have adsmooth_diff3D compiled and then:',
                0165      &      'uncomment the line below and comment the STOP'
                0166            CALL ALL_PROC_DIE( myThid )
                0167            STOP 'ABNORMAL END: S/R smooth_filtervar3D'
                0168 c          CALL adsmooth_diff3D( smoothTmpFld, nbt_in, myThid )
0c3d35c9cd Gael*0169
                0170 c division by sqrt(volume)*sqrt(volume) [1 to end adj, 1 to begin fwd]
9f5240b52a Jean*0171            DO bj=myByLo(myThid),myByHi(myThid)
                0172             DO bi=myBxLo(myThid),myBxHi(myThid)
                0173              DO k=1,Nr
                0174               DO j=1,sNy
                0175                DO i=1,sNx
0c3d35c9cd Gael*0176 c division by ~sqrt(volume):
9f5240b52a Jean*0177                 smoothTmpFld(i,j,k,bi,bj) = smoothTmpFld(i,j,k,bi,bj)
                0178      &                     *(recip_rA(i,j,bi,bj)*recip_drF(k))
                0179                ENDDO
                0180               ENDDO
                0181              ENDDO
                0182             ENDDO
                0183            ENDDO
0c3d35c9cd Gael*0184
                0185 c coming out of adjoint part: overlaps are 0
                0186 c going in fwd part: we need to fill them up
9f5240b52a Jean*0187            CALL EXCH_XYZ_RL( smoothTmpFld , myThid )
0c3d35c9cd Gael*0188
                0189 c fwd:
9f5240b52a Jean*0190            CALL SMOOTH_DIFF3D( smoothTmpFld, nbt_in, myThid )
0c3d35c9cd Gael*0191
                0192 c convert variance to normalization factor:
9f5240b52a Jean*0193            DO bj=myByLo(myThid),myByHi(myThid)
                0194             DO bi=myBxLo(myThid),myBxHi(myThid)
                0195              DO k=1,Nr,dkLoc
                0196               DO j=jj,sNy,djLoc
                0197                DO i=ii,sNx,diLoc
                0198                 IF (smoothTmpFld(i,j,k,bi,bj).NE.0.) THEN
                0199                  smooth3Dnorm(i,j,k,bi,bj) =
                0200      &                    1/SQRT(smoothTmpFld(i,j,k,bi,bj))
                0201                 ENDIF
                0202                ENDDO
                0203               ENDDO
                0204              ENDDO
                0205             ENDDO
                0206            ENDDO
                0207
                0208           ENDDO      !DO ii=1,diLoc
                0209          ENDDO      !DO jj=1,djLoc
                0210         ENDDO      !DO kk=1,dkLoc
                0211
                0212        ELSEIF (smooth3Dfilter(smoothOpNb).EQ.1) THEN
0c3d35c9cd Gael*0213 c compute the normalization matrix using the approximate method
                0214
9f5240b52a Jean*0215         DO bj=myByLo(myThid),myByHi(myThid)
                0216          DO bi=myBxLo(myThid),myBxHi(myThid)
                0217           DO k=1,Nr
                0218            DO j=1-OLy,sNy+OLy
                0219             DO i=1-OLx,sNx+OLx
                0220              smoothTmpMean(i,j,k,bi,bj) = 0. _d 0
                0221              smoothTmpVar(i,j,k,bi,bj)  = 0. _d 0
                0222             ENDDO
                0223            ENDDO
0c3d35c9cd Gael*0224           ENDDO
                0225          ENDDO
                0226         ENDDO
                0227
                0228 c initialize random number generator
9f5240b52a Jean*0229         smoothTmpFld(1,1,1,1,1) = port_rand(1.d0)
                0230         nbRand = 1000
0c3d35c9cd Gael*0231
9f5240b52a Jean*0232         DO ii=1,nbRand
                0233          WRITE(standardMessageUnit,'(A,I4,A,I4)')
                0234      &   'smooth_filtervar3D: ', ii, ' members done out of ', nbRand
0c3d35c9cd Gael*0235
                0236 c fill smoothTmpFld with random numbers:
9f5240b52a Jean*0237          DO bj=myByLo(myThid),myByHi(myThid)
                0238           DO bi=myBxLo(myThid),myBxHi(myThid)
                0239            DO k=1,Nr
                0240             DO j=1-OLy,sNy+OLy
                0241              DO i=1-OLx,sNx+OLx
                0242               smoothTmpFld(i,j,k,bi,bj) = 0. _d 0
                0243               IF (smooth3dmask(i,j,k,bi,bj).NE.0) THEN
                0244                smoothTmpFld(i,j,k,bi,bj)=port_rand_norm()
                0245               ENDIF
0c3d35c9cd Gael*0246 c division by sqrt(volume):
9f5240b52a Jean*0247               smoothTmpFld(i,j,k,bi,bj) = smoothTmpFld(i,j,k,bi,bj)
                0248      &                    *SQRT(recip_rA(i,j,bi,bj)*recip_drF(k))
                0249              ENDDO
                0250             ENDDO
                0251            ENDDO
0c3d35c9cd Gael*0252           ENDDO
                0253          ENDDO
                0254
9f5240b52a Jean*0255          CALL EXCH_XYZ_RL( smoothTmpFld, myThid )
0c3d35c9cd Gael*0256
                0257 c smooth random number field
9f5240b52a Jean*0258          CALL SMOOTH_DIFF3D( smoothTmpFld, nbt_in, myThid )
0c3d35c9cd Gael*0259
                0260 c accumulate statistics (to compute the variance later)
9f5240b52a Jean*0261          DO bj=myByLo(myThid),myByHi(myThid)
                0262           DO bi=myBxLo(myThid),myBxHi(myThid)
                0263            DO k=1,Nr
                0264             DO j=1-OLy,sNy+OLy
                0265              DO i=1-OLx,sNx+OLx
                0266               smoothTmpVar(i,j,k,bi,bj) = smoothTmpVar(i,j,k,bi,bj)
                0267      &         + smoothTmpFld(i,j,k,bi,bj)
                0268      &          *smoothTmpFld(i,j,k,bi,bj)/nbRand
                0269               smoothTmpMean(i,j,k,bi,bj) = smoothTmpMean(i,j,k,bi,bj)
                0270      &         + smoothTmpFld(i,j,k,bi,bj)/nbRand
                0271              ENDDO
                0272             ENDDO
                0273            ENDDO
0c3d35c9cd Gael*0274           ENDDO
                0275          ENDDO
                0276
9f5240b52a Jean*0277 C-      ii loop end
                0278         ENDDO
0c3d35c9cd Gael*0279
                0280 c compute variance and convert it to normalization factor:
9f5240b52a Jean*0281         DO bj=myByLo(myThid),myByHi(myThid)
                0282          DO bi=myBxLo(myThid),myBxHi(myThid)
                0283           DO k=1,Nr
                0284            DO j=1-OLy,sNy+OLy
                0285             DO i=1-OLx,sNx+OLx
                0286              IF (smooth3dmask(i,j,k,bi,bj).NE.0) THEN
                0287               smooth3Dnorm(i,j,k,bi,bj) = 1/SQRT(
                0288      &           nbRand/(nbRand-1)
                0289      &          *( smoothTmpVar(i,j,k,bi,bj)
                0290      &           - smoothTmpMean(i,j,k,bi,bj)*smoothTmpMean(i,j,k,bi,bj)
                0291      &           )                              )
                0292              ENDIF
                0293             ENDDO
                0294            ENDDO
0c3d35c9cd Gael*0295           ENDDO
                0296          ENDDO
                0297         ENDDO
                0298
9f5240b52a Jean*0299 C      end smooth3Dfilter() if block
                0300        ENDIF
0c3d35c9cd Gael*0301
                0302 c write smooth3Dnorm_3D to file:
de57a2ec4b Mart*0303        WRITE(fnamegeneric,'(2A,I3.3)')
f9d7cbfb72 Ou W*0304      &       smoothDir(1:IL),
9f5240b52a Jean*0305      &       'smooth3Dnorm', smoothOpNb
                0306        CALL WRITE_REC_3D_RL( fnamegeneric, smoothprec,
                0307      &                       Nr, smooth3Dnorm, 1, 1, myThid )
                0308        CALL EXCH_XYZ_RL( smooth3Dnorm,  myThid )
0c3d35c9cd Gael*0309
9f5240b52a Jean*0310 C     end smooth3Dfilter() <> 0 if block
0c3d35c9cd Gael*0311       ENDIF
                0312
                0313       RETURN
                0314       END