pkg/smooth/smooth_filtervar2d.F

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