pkg/ebm/ebm_wind_perturb.F

d8206d87ee Patr*0001 #include "EBM_OPTIONS.h"
                0002
                0003 CStartOfInterface
b08554040b Patr*0004       SUBROUTINE EBM_WIND_PERTURB( myTime, myIter, myThid )
d7d2f03550 Jean*0005 C     *==========================================================*
                0006 C     | S/R EBM_WIND_PERTURB
                0007 C     | o Calculated random wind perturbations.
                0008 C     *==========================================================*
d8206d87ee Patr*0009       IMPLICIT NONE
                0010
                0011 C     == Global data ==
                0012 #include "SIZE.h"
                0013 #include "EEPARAMS.h"
                0014 #include "PARAMS.h"
                0015 #include "GRID.h"
                0016 #include "DYNVARS.h"
                0017 #include "FFIELDS.h"
                0018 #ifdef ALLOW_EBM
                0019 # include "EBM.h"
                0020 #endif
                0021
                0022 C     == Routine arguments ==
b08554040b Patr*0023       _RL    myTime
                0024       INTEGER myIter
                0025       INTEGER myThid
d8206d87ee Patr*0026 CEndOfInterface
                0027
                0028 #ifdef ALLOW_EBM
                0029 # ifdef EBM_WIND_PERT
                0030
                0031 C     == Local variables ==
                0032 C     Loop counters
b08554040b Patr*0033       INTEGER i, j, bi, bj
d7d2f03550 Jean*0034       _RS ya(1-OLy:sNy+OLy), ya2(1-OLy:sNy+OLy)
                0035       _RS xa(1-OLx:sNx+OLx), xa2(1-OLx:sNx+OLx)
d8206d87ee Patr*0036       _RS y(1-OLy:sNy+OLy), x(1-OLx:sNx+OLx)
                0037       _RS temp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0038       _RS temp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0039       _RS stdev(1-OLy:sNy+OLy)
                0040       _RS std(1:40)
                0041       data std /0.030, 0.035, 0.045, 0.053, 0.059, 0.060, 0.056,
                0042      &          0.048, 0.041, 0.038, 0.034, 0.029, 0.023, 0.018,
                0043      &          0.016, 0.015, 0.013, 0.011, 0.008, 0.005, 0.005,
                0044      &          0.005, 0.008, 0.011, 0.014, 0.014, 0.017, 0.019,
                0045      &          0.023, 0.029, 0.032, 0.038, 0.048, 0.058, 0.065,
                0046      &          0.067, 0.063, 0.060, 0.062, 0.064 /
                0047
                0048
b08554040b Patr*0049       DO bj=myByLo(myThid),myByHi(myThid)
                0050        DO bi=myBxLo(myThid),myBxHi(myThid)
                0051
d8206d87ee Patr*0052       DO j = 1-OLy, sNy+OLy
                0053          y(j) = 0.0
                0054          ya(j) = 0.0
                0055          ya2(j) =  0.0
                0056          stdev(j) = 0.0
                0057       ENDDO
                0058       DO i = 1-OLx, sNx+OLx
                0059          x(i) = 0.0
                0060          xa(i) = 0.0
                0061          xa2(i) =  0.0
                0062       ENDDO
                0063       DO i = 1-OLx, sNx+OLx
                0064         DO j = 1-OLy, sNy+OLy
                0065          temp(i,j) = 0.0
                0066          temp2(i,j) = 0.0
                0067         ENDDO
                0068       ENDDO
                0069       DO j = 1, sNy
                0070          stdev(j) = std(j)
                0071       ENDDO
                0072
                0073 cph   Generate random numbers
                0074 cph   Need to get this from somewhere!
                0075       call random_number (temp)
                0076
                0077 C     interpolation in first dimension
                0078 C     scaling to get a process with a standard deviation of 1
                0079       DO j = jMin, jMax
                0080        DO i = iMin, iMax
                0081          temp(i,j) = 1.73*(2.0*temp(i,j) - 1.0)
                0082        ENDDO
                0083       ENDDO
                0084
                0085       DO j = jMin, jMax
                0086        DO i = iMin, iMax
                0087          x(i) = i
                0088          xa(i) = x(i) - MOD(x(i),10.0)
                0089          xa2(i) = xa(i)+10.0
                0090         if ( xa2(i) .gt. sNx+Olx ) then
                0091            xa2(i) = 0.0
                0092          endif
                0093          temp2(i,j) = 0.1*( (x(i)-xa(i))*temp(INT(xa2(i)),j)+
                0094      &        (10.0-x(i)+xa(i))*temp(INT(xa(i)),j) )
                0095        ENDDO
                0096       ENDDO
                0097
                0098 C     interpolation in second dimension
                0099 C     multiplication with observation zonal wind stress standard deviation
                0100 C     add AR1 process
                0101       DO i = iMin, iMax
                0102        DO j = jMin, jMax
                0103          y(j) = j
                0104          ya(j) = y(j) - MOD(y(j),6.0)
                0105          ya2(j) = ya(j)+6.0
                0106          if ( ya2(j) .gt. sNy+Oly ) then
                0107             ya2(j) = 0.0
                0108          endif
d7d2f03550 Jean*0109 c     time lag correlation coefficient, use 0.75 for temperature timescale,
d8206d87ee Patr*0110 c     0.98 for the momentum timescale.
                0111          winPert(i,j,bi,bj) = maskW(i,j,k,bi,bj)*
d7d2f03550 Jean*0112      &        (1.0/1.66)*(0.75*winPert(i,j,bi,bj) +
d8206d87ee Patr*0113      &        stdev(j)*(1.0/6.0)*
                0114      &        ((y(j)-ya(j))*temp2(i,INT(ya2(j)))+
                0115      &        (6.0-y(j)+ya(j))*temp2(i,INT(ya(j)))))
                0116        ENDDO
                0117       ENDDO
                0118
b08554040b Patr*0119        ENDDO
                0120       ENDDO
                0121
d8206d87ee Patr*0122 C      CALL PLOT_FIELD_XYRS( winPert, 'winPert',1,myThid)
                0123
6637358eea Jean*0124       _EXCH_XY_RS(winPert  , myThid )
d8206d87ee Patr*0125
d7d2f03550 Jean*0126       DO bj=myByLo(myThid),myByHi(myThid)
                0127        DO bi=myBxLo(myThid),myBxHi(myThid)
                0128          DO j = 1-OLy, sNy+OLy
                0129           DO i = 1-OLx, sNx+OLx
                0130             fu(i,j,bi,bj) = fu(i,j,bi,bj)
                0131      &                    + winPert(i,j,bi,bj)*rUnit2mass
                0132      &                     *drF(1)*hFacW(i,j,1,bi,bj)
                0133           ENDDO
                0134          ENDDO
                0135        ENDDO
                0136       ENDDO
                0137
d8206d87ee Patr*0138 # endif /* EBM_WIND_PERT */
                0139 #endif /* ALLOW_EBM */
                0140
d7d2f03550 Jean*0141       RETURN
                0142       END