/[MITgcm]/MITgcm/pkg/kl10/kl10_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/kl10/kl10_calc.F

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

revision 1.1 by jmc, Wed Jul 30 03:28:05 2014 UTC revision 1.2 by jmc, Wed Jul 30 19:52:38 2014 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: KL10_CALC  C !ROUTINE: KL10_CALC
8    
9  C !INTERFACE: =======================================================  C !INTERFACE: =======================================================
10         subroutine KL10_CALC(        SUBROUTINE KL10_CALC(
11       I     bi, bj, myTime, myThid )       I                bi, bj, sigmaR, myTime, myIter, myThid )
12    
13  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     /==========================================================\  C     *==========================================================*
15  C     | SUBROUTINE KL10_CALC                                     |  C     | SUBROUTINE KL10_CALC                                     |
16  C     | o Compute all KL10 fields defined in KL10.h              |  C     | o Compute all KL10 fields defined in KL10.h              |
17  C     |==========================================================|  C     *==========================================================*
18  C     | This subroutine is based on SPEM code                    |  C     | This subroutine is based on SPEM code                    |
19  C     \==========================================================/  C     *==========================================================*
20        IMPLICIT NONE  C \ev
21  C  
22  C--------------------------------------------------------------------  C--------------------------------------------------------------------
23    
24  C JMK  C JMK
25  C global parameters updated by kl_calc  C global parameters updated by kl_calc
26  C     KLviscAz   - KL eddy viscosity coefficient              (m^2/s)  C     KLviscAz  :: KL eddy viscosity coefficient              (m^2/s)
27  C     KLdiffKzT  - KL diffusion coefficient for temperature   (m^2/s)  C     KLdiffKzT :: KL diffusion coefficient for temperature   (m^2/s)
 C  
 C \ev  
28    
29  C !USES: ============================================================  C !USES: ============================================================
30          IMPLICIT NONE
31  #include "SIZE.h"  #include "SIZE.h"
32  #include "EEPARAMS.h"  #include "EEPARAMS.h"
33  #include "PARAMS.h"  #include "PARAMS.h"
# Line 46  C !USES: =============================== Line 45  C !USES: ===============================
45    
46  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
47  c Routine arguments  c Routine arguments
48  c     bi, bj - array indices on which to apply calculations  C     bi, bj :: Current tile indices
49  c     myTime - Current time in simulation  C     sigmaR :: Vertical gradient of iso-neutral density
50    C     myTime :: Current time in simulation
51         INTEGER bi, bj  C     myIter :: Current time-step number
52         INTEGER myThid  C     myThid :: My Thread Id number
53         _RL     myTime        INTEGER bi, bj
54          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
55          _RL     myTime
56          INTEGER myIter
57          INTEGER myThid
58    
59  #ifdef ALLOW_KL10  #ifdef ALLOW_KL10
   
60  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
61  c Local constants  c Local constants
62  C     imin, imax, jmin, jmax  - array computation indices  C     iMin, iMax, jMin, jMax :: array computation indices
63  C     RiNumber                - Richardson Number  C     RiNumber               :: Richardson Number
64        INTEGER I, J, K, Km1, JJ        INTEGER I, J, K, Km1, JJ
65        INTEGER   iMin ,iMax ,jMin ,jMax,di        INTEGER   iMin ,iMax ,jMin ,jMax,di
66        _RL     denom, KLviscTmp, buoyFreq,rhot,tempu,tempv,oldK,Ri        _RL     denom, KLviscTmp, buoyFreq,rhot,tempu,tempv,oldK,Ri
67        _RL     b0, buoyFreqf, buoyFreqc, KLviscold,zsum,zsums        _RL     b0, buoyFreqf, buoyFreqc, KLviscold,zsum,zsums
68        _RL     rhoS(1:Nr),RS(1:Nr)        _RL     rhoS(1:Nr),RS(1:Nr)
69        _RL     dzc,dzp,ec,ep,es,epss(-1:0),epsw(-1:0),dz,KTemp        _RL     dzc,dzp,ec,ep,es,epss(-1:0),epsw(-1:0),dz,KTemp
70  !      _RL     bF(1:Nr)  c     _RL     bF(1:Nr)
71  !      _RL     theta_mcb(1:Nr),theta_mcb3(1:Nr)  c     _RL     theta_mcb(1:Nr),theta_mcb3(1:Nr)
72  C     === Local variables ===  C     === Local variables ===
73  C     msgBuf      - Informational/error message buffer  C     msgBuf     :: Informational/error message buffer
74        CHARACTER*(1024) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
   
75  CEOP  CEOP
76    
77        iMin = 2-OLx        iMin = 2-OLx
78        iMax = sNx+OLx-1        iMax = sNx+OLx-1
79        jMin = 2-OLy        jMin = 2-OLy
# Line 86  CEOP Line 88  CEOP
88              RS(1)=rC(1)              RS(1)=rC(1)
89    
90              KLeps(I-1,J-1,1,bi,bj)=0.0              KLeps(I-1,J-1,1,bi,bj)=0.0
91  C           eps(k-1) = (dz(k-1)*eps0(k-1) +dz(k)*eps0(k))/(dz(k-1)+dz(k))  c           eps(k-1) = (dz(k-1)*eps0(k-1) +dz(k)*eps0(k))/(dz(k-1)+dz(k))
92              ep = 0.0              ep = 0.0
93              dzp = 0.0              dzp = 0.0
94    
# Line 99  C the density came from, RS. Line 101  C the density came from, RS.
101       $              ,bj),totPhiHyd(I,J,K,bi,bj),rhot,myThid )       $              ,bj),totPhiHyd(I,J,K,bi,bj),rhot,myThid )
102                 rhoS(K)=rhot                 rhoS(K)=rhot
103                 RS(K)=rC(K)                 RS(K)=rC(K)
104  C$$$               WRITE(msgBuf, '(A,I10.10,A,E10.4,A,E10.4)') 'Hellok ', K  c$$$               WRITE(msgBuf, '(A,I10.10,A,E10.4,A,E10.4)') 'Hellok ', K
105  C$$$     $              -1,' ',theta(I,J,K,bi,bj),' ',rhot  c$$$     $              -1,' ',theta(I,J,K,bi,bj),' ',rhot
106  C$$$               CALL PRINT_MESSAGE(msgBuf, standardMessageUnit,  c$$$               CALL PRINT_MESSAGE(msgBuf, standardMessageUnit,
107  C$$$     &              SQUEEZE_RIGHT , 1)  c$$$     &              SQUEEZE_RIGHT , 1)
108                 IF ( (rhoS(K).LT.rhoS(K-1)).AND.(maskC(I,J,K,bi                 IF ( (rhoS(K).LT.rhoS(K-1)).AND.(maskC(I,J,K,bi
109       &              ,bj).GT.0)) THEN       &              ,bj).GT.0)) THEN
110                    JJ=K-1                    JJ=K-1
111                    DO WHILE ( (JJ.GT.0).AND.(rhoS(K).LT.rhoS(JJ)) )                    DO WHILE ( (JJ.GT.0).AND.(rhoS(K).LT.rhoS(JJ)) )
112  C                     write(*,*) K,JJ,rhoS(K),rhoS(JJ)  c                    write(*,*) K,JJ,rhoS(K),rhoS(JJ)
113                       JJ=JJ-1                       JJ=JJ-1
114                    ENDDO                    ENDDO
115                    rhoS(JJ+1:K)=cshift(rhoS(JJ+1:K),-1)                    rhoS(JJ+1:K)=cshift(rhoS(JJ+1:K),-1)
# Line 125  C diffKrNrT(K) = viscArNr(K) = backgroun Line 127  C diffKrNrT(K) = viscArNr(K) = backgroun
127  C N at surface = zero or uses gradient  C N at surface = zero or uses gradient
128              b0 = MAX(-gravity*mass2rUnit*              b0 = MAX(-gravity*mass2rUnit*
129       &              (rhoS(1) - rhoS(2))*recip_drC(2),0. _d 0)       &              (rhoS(1) - rhoS(2))*recip_drC(2),0. _d 0)
130  C            b0 = 0.  c           b0 = 0.
131              DO di=-1,0              DO di=-1,0
132                 epss(di)=0.0                 epss(di)=0.0
133                 epsw(di)=0.0                 epsw(di)=0.0
# Line 228  C add the u and v dissipations: Line 230  C add the u and v dissipations:
230  C Note this ec is defined on cell faces K=2..NR at the center of the  C Note this ec is defined on cell faces K=2..NR at the center of the
231  C cells (i.e. at XC), so its above the density variables.  C cells (i.e. at XC), so its above the density variables.
232  C  C
233  C So to get at the center of the cells, just average this one and the previous one.  And its a true average because the  C So to get at the center of the cells, just average this one and the previous one.
234    C  And its a true average because the
235    
236                    KLeps(I-1,J-1,Km1,bi,bj) = 0.5*(ep+ec)                    KLeps(I-1,J-1,Km1,bi,bj) = 0.5*(ep+ec)
237                    IF (Km1.EQ.1) THEN                    IF (Km1.EQ.1) THEN
# Line 258  C           ENDDO I Line 261  C           ENDDO I
261    
262        RETURN        RETURN
263        END        END
   
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22