/[MITgcm]/MITgcm/verification/hs94.cs-32x32x5/code/ini_theta.F
ViewVC logotype

Diff of /MITgcm/verification/hs94.cs-32x32x5/code/ini_theta.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by adcroft, Wed Jun 6 19:46:43 2001 UTC revision 1.4 by jmc, Fri Jul 6 22:13:37 2001 UTC
# Line 42  C     I,J,K Line 42  C     I,J,K
42        INTEGER bi, bj        INTEGER bi, bj
43        INTEGER I, J, K        INTEGER I, J, K
44        _RL     term1,term2,thetaLim,thetaEq        _RL     term1,term2,thetaLim,thetaEq
       _RL     thKappa  
45    
46        _BARRIER        _BARRIER
47    
# Line 50  C     I,J,K Line 49  C     I,J,K
49  c     CALL SRAND( J )  c     CALL SRAND( J )
50    
51        IF ( hydrogThetaFile .EQ. ' ' ) THEN        IF ( hydrogThetaFile .EQ. ' ' ) THEN
52  C--    Initialise temperature field to Held & Saurez equilibrium theta  C--    Initialise temperature field to Held & Suarez equilibrium theta
53         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
54          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
55           DO K=1,Nr           DO K=1,Nr
56            Ro_SeaLevel=1. _d 5            thetaLim = 200. _d 0/((rC(K)/atm_po)**atm_kappa)
           thKappa = 2. _d 0/7. _d 0  
           thetaLim = 200. _d 0/((rC(K)/Ro_SeaLevel)**thKappa)  
57            DO J=1,sNy            DO J=1,sNy
58             DO I=1,sNx             DO I=1,sNx
59             term1=60. _d 0*(sin(yC(I,J,bi,bj)*deg2rad)**2)             term1=60. _d 0*(sin(yC(I,J,bi,bj)*deg2rad)**2)
60             term2=10. _d 0*log((rC(K)/Ro_SeaLevel))             term2=10. _d 0*log((rC(K)/atm_po))
61       &              *(cos(yC(I,J,bi,bj)*deg2rad)**2)       &              *(cos(yC(I,J,bi,bj)*deg2rad)**2)
62             thetaEq=315. _d 0-term1-term2             thetaEq=315. _d 0-term1-term2
63              theta(I,J,K,bi,bj) = MAX( thetaLim, thetaEq )              theta(I,J,K,bi,bj) = MAX( thetaLim, thetaEq )
# Line 72  c           theta(I,J,K,bi,bj) = tRef(K) Line 69  c           theta(I,J,K,bi,bj) = tRef(K)
69           ENDDO           ENDDO
70          ENDDO          ENDDO
71         ENDDO         ENDDO
72    #ifdef ALLOW_ZONAL_FILT
73    C--   Zonal FFT filter initial conditions
74         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
75          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
76           DO K=1,Nr           DO K=1,Nr
77            DO J=1,sNy            DO J=1,sNy
 #ifdef ALLOW_ZONAL_FILT  
 C--   Zonal FFT filter initial conditions  
78             CALL ZONAL_FILTER(             CALL ZONAL_FILTER(
79       U      theta, hFacC,       U      theta, hFacC,
80       I      1, sNy, k, k, bi, bj, 1, myThid)       I      1, sNy, k, k, bi, bj, 1, myThid)
 #endif /* INCLUDE_LAT_CIRC_FFT_FILTER_CODE */  
81            ENDDO            ENDDO
82           ENDDO           ENDDO
83          ENDDO          ENDDO
84         ENDDO         ENDDO
85    #endif /* INCLUDE_LAT_CIRC_FFT_FILTER_CODE */
86        ELSE        ELSE
87         _BEGIN_MASTER( myThid )         _BEGIN_MASTER( myThid )
88         CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )         CALL READ_FLD_XYZ_RL( hydrogThetaFile, ' ', theta, 0, myThid )
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22