/[MITgcm]/MITgcm/pkg/layers/layers_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/layers/layers_calc.F

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

revision 1.8 by dfer, Thu May 12 15:09:54 2011 UTC revision 1.15 by rpa, Wed Oct 17 18:49:15 2012 UTC
# Line 16  C !INTERFACE: Line 16  C !INTERFACE:
16  C !DESCRIPTION:  C !DESCRIPTION:
17  C ===================================================================  C ===================================================================
18  C     Calculate the transport in isopycnal layers.  C     Calculate the transport in isopycnal layers.
19  C     This is the meat of the LAYERS package.  C     This was the meat of the LAYERS package, which
20    C     has been moved to S/R LAYERS_FLUXCALC.F
21  C ===================================================================  C ===================================================================
22    
23  C !USES:  C !USES:
# Line 46  CEOP Line 47  CEOP
47  C !LOCAL VARIABLES:  C !LOCAL VARIABLES:
48  C     bi, bj   :: tile indices  C     bi, bj   :: tile indices
49  C     i,j      :: horizontal indices  C     i,j      :: horizontal indices
50    C     iLa      :: layer coordinate index
51  C     k        :: vertical index for model grid  C     k        :: vertical index for model grid
52  C     kci      :: index from CellIndex  
53  C     kg       :: index for looping though layers_G        INTEGER bi, bj, iLa
54  C     kk       :: vertical index for ZZ (fine) grid  #ifdef LAYERS_PRHO_REF
55  C     kgu,kgv  :: vertical index for isopycnal grid        INTEGER i, j, k
 C     TatV     :: temperature at U point  
 C     TatV     :: temperature at V point  
   
       INTEGER bi, bj  
       INTEGER i,j,k,kk,kg,kci  
       INTEGER kgu(sNx+1,sNy+1), kgv(sNx+1,sNy+1)  
       _RL TatU, TatV  
       CHARACTER*(MAX_LEN_MBUF) msgBuf  
 #if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)  
       INTEGER kcip1  
       _RL delPsi, maskp1  
56  #endif  #endif
57    
58  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
59    
60  C --- The tile loops        DO iLa=1,layers_maxNum
61        DO bj=myByLo(myThid),myByHi(myThid)        
62        DO bi=myBxLo(myThid),myBxHi(myThid)        IF (layers_num(iLa) .EQ. 1) THEN
63           CALL LAYERS_FLUXCALC( uVel,vVel,theta,iLa,
64  C     Initialize the search indices       &  layers_UH(1-OLx,1-OLy,1,1,1,iLa),
65        DO j = 1,sNy+1       &  layers_VH(1-OLx,1-OLy,1,1,1,iLa),
66          DO i = 1,sNx+1       &  layers_Hw(1-OLx,1-OLy,1,1,1,iLa),
67  C       The temperature index (layer_G) goes from cold to warm.       &  layers_Hs(1-OLx,1-OLy,1,1,1,iLa),
68  C       The water column goes from warm (k=1) to cold (k=Nr).       &  layers_PIw(1-OLx,1-OLy,1,1,1,iLa),
69  C       So initialize the search with the warmest value.       &  layers_PIs(1-OLx,1-OLy,1,1,1,iLa),
70            kgu(i,j) = Nlayers       &  layers_U(1-OLx,1-OLy,1,1,1,iLa),
71            kgv(i,j) = Nlayers       &  layers_V(1-OLx,1-OLy,1,1,1,iLa),
72          ENDDO       &  myThid)
73        ENDDO        ELSEIF (layers_num(iLa) .EQ. 2) THEN
74           CALL LAYERS_FLUXCALC( uVel,vVel,salt,iLa,
75  C     Reset the arrays       &  layers_UH(1-OLx,1-OLy,1,1,1,iLa),
76        DO kg=1,Nlayers       &  layers_VH(1-OLx,1-OLy,1,1,1,iLa),
77         DO j = 1,sNy+1       &  layers_Hw(1-OLx,1-OLy,1,1,1,iLa),
78          DO i = 1,sNx+1       &  layers_Hs(1-OLx,1-OLy,1,1,1,iLa),
79  #ifdef LAYERS_UFLUX       &  layers_PIw(1-OLx,1-OLy,1,1,1,iLa),
80           layers_UFlux(i,j,kg,bi,bj) = 0. _d 0       &  layers_PIs(1-OLx,1-OLy,1,1,1,iLa),
81  #ifdef LAYERS_THICKNESS       &  layers_U(1-OLx,1-OLy,1,1,1,iLa),
82           layers_HU(i,j,kg,bi,bj) = 0. _d 0       &  layers_V(1-OLx,1-OLy,1,1,1,iLa),
83  #endif /* LAYERS_THICKNESS */       &  myThid)
84  #endif /* LAYERS_UFLUX */        ELSEIF (layers_num(iLa) .EQ. 3) THEN
85  #ifdef LAYERS_VFLUX  #ifdef LAYERS_PRHO_REF
86           layers_VFlux(i,j,kg,bi,bj) = 0. _d 0  C     For layers_num(iLa) = 3, calculate the potential density referenced to
87  #ifdef LAYERS_THICKNESS  C     the model level given by layers_krho.
88           layers_HV(i,j,kg,bi,bj) = 0. _d 0         DO bj=myByLo(myThid),myByHi(myThid)
89  #endif /* LAYERS_THICKNESS */         DO bi=myBxLo(myThid),myBxHi(myThid)
90  #endif /* LAYERS_VFLUX */         DO k = 1,Nr
91            CALL FIND_RHO_2D( 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,
92         &     layers_krho(iLa),
93         &     theta(1-OLx,1-OLy,k,bi,bj),
94         &     salt(1-OLx,1-OLy,k,bi,bj),
95         &     prho(1-OLx,1-OLy,k,bi,bj,iLa),
96         &     k, bi, bj, myThid )
97            DO j = 1-OLy,sNy+OLy
98             DO i = 1-OLx,sNx+OLx
99               prho(i,j,k,bi,bj,iLa) = rhoConst + prho(i,j,k,bi,bj,iLa)
100             ENDDO
101          ENDDO          ENDDO
102         ENDDO         ENDDO
103        ENDDO         ENDDO
104           ENDDO
105  C      _RL  theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)         CALL LAYERS_FLUXCALC( uVel,vVel,
106  C     Sometimes it is done this way       &  prho(1-OLx,1-OLy,1,1,1,iLa),iLa,
107  C      DO j=1-Oly+1,sNy+Oly-1       &  layers_UH(1-OLx,1-OLy,1,1,1,iLa),
108  C       DO i=1-Olx+1,sNx+Olx-1       &  layers_VH(1-OLx,1-OLy,1,1,1,iLa),
109        DO kk=1,NZZ       &  layers_Hw(1-OLx,1-OLy,1,1,1,iLa),
110         k = MapIndex(kk)       &  layers_Hs(1-OLx,1-OLy,1,1,1,iLa),
111         kci = CellIndex(kk)       &  layers_PIw(1-OLx,1-OLy,1,1,1,iLa),
112         DO j = 1,sNy+1       &  layers_PIs(1-OLx,1-OLy,1,1,1,iLa),
113          DO i = 1,sNx+1       &  layers_U(1-OLx,1-OLy,1,1,1,iLa),
114         &  layers_V(1-OLx,1-OLy,1,1,1,iLa),    
115  #ifdef LAYERS_UFLUX       &  myThid)
 C ------ Find theta at the U point (west) on the fine Z grid  
          IF (LAYER_nb .EQ. 1) THEN  
          TatU = MapFact(kk) *  
      &    0.5 _d 0 * (theta(i-1,j,k,bi,bj)+theta(i,j,k,bi,bj)) +  
      &    (1-MapFact(kk)) *  
      &    0.5 _d 0 * (theta(i-1,j,k+1,bi,bj)+theta(i,j,k+1,bi,bj))  
          ELSEIF (LAYER_nb .EQ. 2) THEN  
          TatU = MapFact(kk) *  
      &    0.5 _d 0 * (salt(i-1,j,k,bi,bj)+salt(i,j,k,bi,bj)) +  
      &    (1-MapFact(kk)) *  
      &    0.5 _d 0 * (salt(i-1,j,k+1,bi,bj)+salt(i,j,k+1,bi,bj))  
          ENDIF  
   
 C ------ Now that we know T everywhere, determine the binning.  
   
          IF (TatU .GE. layers_G(Nlayers)) THEN  
 C        the point is in the hottest bin or hotter  
           kgu(i,j) = Nlayers  
          ELSE IF (TatU .LT. layers_G(2)) THEN  
 C        the point is in the coldest bin or colder  
           kgu(i,j) = 1  
          ELSE IF ( (TatU .GE. layers_G(kgu(i,j)))  
      &    .AND. (TatU .LT. layers_G(kgu(i,j)+1)) ) THEN  
 C        already on the right bin -- do nothing  
          ELSE IF (TatU .GE. layers_G(kgu(i,j))) THEN  
 C        have to hunt for the right bin by getting hotter  
           DO WHILE (TatU .GE. layers_G(kgu(i,j)+1))  
            kgu(i,j) = kgu(i,j) + 1  
           ENDDO  
 C         now layers_G(kgu(i,j)+1) < TatU <= layers_G(kgu(i,j)+1)  
          ELSE IF (TatU .LT. layers_G(kgu(i,j)+1)) THEN  
 C        have to hunt for the right bin by getting colder  
           DO WHILE (TatU .LT. layers_G(kgu(i,j)))  
            kgu(i,j) = kgu(i,j) - 1  
           ENDDO  
 C         now layers_G(kgu(i,j)+1) <= TatU < layers_G(kgu(i,j)+1)  
          ELSE  
 C        that should have covered all the options  
           WRITE(msgBuf,'(A,1E14.6)')  
      &     'S/R LAYERS_CALC: Couldnt find a bin in layers_G for TatU=',  
      &     TatU  
           CALL PRINT_ERROR( msgBuf, myThid )  
           STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'  
          END IF  
   
 C ------ Augment the bin values  
          layers_UFlux(i,j,kgu(i,j),bi,bj) =  
      &    layers_UFlux(i,j,kgu(i,j),bi,bj) +  
      &    dZZf(kk) * uVel(i,j,kci,bi,bj) * hFacW(i,j,kci,bi,bj)  
   
 #if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)  
          IF ( GM_AdvForm .AND. useBOLUS ) THEN  
            kcip1 = MIN(kci+1,Nr)  
            maskp1 = 1.  
            IF (kci.GE.Nr) maskp1 = 0.  
            delPsi = GM_PsiX(i,j,kcip1,bi,bj)*maskp1  
      &            - GM_PsiX(i,j, kci, bi,bj)  
            layers_UFlux(i,j,kgu(i,j),bi,bj) =  
      &      layers_UFlux(i,j,kgu(i,j),bi,bj)  
      &      + delPsi*recip_drF(kci)*_recip_hFacW(i,j,kci,bi,bj)  
      &      * dZZf(kk)*hFacW(i,j,kci,bi,bj)  
          ENDIF  
 #endif  
   
 #ifdef LAYERS_THICKNESS  
          layers_HU(i,j,kgu(i,j),bi,bj) = layers_HU(i,j,kgu(i,j),bi,bj)  
      &    + dZZf(kk) * hFacW(i,j,kci,bi,bj)  
 #endif /* LAYERS_THICKNESS */  
   
 #endif /* LAYERS_UFLUX */  
   
 #ifdef LAYERS_VFLUX  
 C ------ Find theta at the V point (south) on the fine Z grid  
          IF (LAYER_nb .EQ. 1) THEN  
          TatV = MapFact(kk) *  
      &    0.5 _d 0 * (theta(i,j-1,k,bi,bj)+theta(i,j,k,bi,bj)) +  
      &    (1-MapFact(kk)) *  
      &    0.5 _d 0 * (theta(i,j-1,k+1,bi,bj)+theta(i,j,k+1,bi,bj))  
          ELSEIF (LAYER_nb .EQ. 2) THEN  
          TatV = MapFact(kk) *  
      &    0.5 _d 0 * (salt(i,j-1,k,bi,bj)+salt(i,j,k,bi,bj)) +  
      &    (1-MapFact(kk)) *  
      &    0.5 _d 0 * (salt(i,j-1,k+1,bi,bj)+salt(i,j,k+1,bi,bj))  
          ENDIF  
   
 C ------ Now that we know T everywhere, determine the binning  
          IF (TatV .GE. layers_G(Nlayers)) THEN  
 C         the point is in the hottest bin or hotter  
           kgv(i,j) = Nlayers  
          ELSE IF (TatV .LT. layers_G(2)) THEN  
 C         the point is in the coldest bin or colder  
           kgv(i,j) = 1  
          ELSE IF ( (TatV .GE. layers_G(kgv(i,j)))  
      &    .AND. (TatV .LT. layers_G(kgv(i,j)+1)) ) THEN  
 C         already on the right bin -- do nothing  
          ELSE IF (TatV .GE. layers_G(kgv(i,j))) THEN  
 C         have to hunt for the right bin by getting hotter  
           DO WHILE (TatV .GE. layers_G(kgv(i,j)+1))  
            kgv(i,j) = kgv(i,j) + 1  
           ENDDO  
 C         now layers_G(kgv(i,j)+1) < TatV <= layers_G(kgv(i,j)+1)  
          ELSE IF (TatV .LT. layers_G(kgv(i,j)+1)) THEN  
 C         have to hunt for the right bin by getting colder  
           DO WHILE (TatV .LT. layers_G(kgv(i,j)))  
            kgv(i,j) = kgv(i,j) - 1  
           ENDDO  
 C         now layers_G(kgv(i,j)+1) <= TatV < layers_G(kgv(i,j)+1)  
          ELSE  
 C         that should have covered all the options  
           WRITE(msgBuf,'(A,1E14.6)')  
      &     'S/R LAYERS_CALC: Couldnt find a bin in layers_G for TatV=',  
      &     TatV  
           CALL PRINT_ERROR( msgBuf, myThid )  
           STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'  
          END IF  
   
 C ------ Augment the bin values  
          layers_VFlux(i,j,kgv(i,j),bi,bj) =  
      &    layers_VFlux(i,j,kgv(i,j),bi,bj)  
      &    + dZZf(kk) * vVel(i,j,kci,bi,bj) * hFacS(i,j,kci,bi,bj)  
   
 #if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)  
          IF ( GM_AdvForm .AND. useBOLUS ) THEN  
            kcip1 = MIN(kci+1,Nr)  
            maskp1 = 1.  
            IF (kci.GE.Nr) maskp1 = 0.  
            delPsi = GM_PsiY(i,j,kcip1,bi,bj)*maskp1  
      &            - GM_PsiY(i,j, kci, bi,bj)  
            layers_VFlux(i,j,kgv(i,j),bi,bj) =  
      &      layers_VFlux(i,j,kgv(i,j),bi,bj)  
      &      + delPsi*recip_drF(kci)*_recip_hFacS(i,j,kci,bi,bj)  
      &      * dZZf(kk)*hFacS(i,j,kci,bi,bj)  
          ENDIF  
116  #endif  #endif
117          ENDIF
118    
119  #ifdef LAYERS_THICKNESS  #ifdef ALLOW_LAYERS_OUTPUT
          layers_HV(i,j,kgv(i,j),bi,bj) = layers_HV(i,j,kgv(i,j),bi,bj)  
      &    + dZZf(kk) * hFacS(i,j,kci,bi,bj)  
 #endif /* LAYERS_THICKNESS */  
   
 #endif /* LAYERS_VFLUX */  
   
         ENDDO  
        ENDDO  
       ENDDO  
120    
121  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
122  C--   Time-average  C--   Time-average
123    cgf layers_maxNum loop and dimension would be needed for
124    cgf the following and tave output to work beyond iLa.EQ.1
125          IF ( iLa.EQ.1 ) THEN
126        IF ( layers_taveFreq.GT.0. ) THEN        IF ( layers_taveFreq.GT.0. ) THEN
127    C --- The tile loops
128          DO bj=myByLo(myThid),myByHi(myThid)
129          DO bi=myBxLo(myThid),myBxHi(myThid)
130    
131  #ifdef LAYERS_UFLUX  #ifdef LAYERS_UFLUX
132           CALL TIMEAVE_CUMULATE( layers_UFlux_T, layers_UFlux, Nlayers,           CALL TIMEAVE_CUMULATE( layers_UH_T, layers_UFlux, Nlayers,
133       &                          deltaTclock, bi, bj, myThid )       &                          deltaTclock, bi, bj, myThid )
134  #ifdef LAYERS_THICKNESS  #ifdef LAYERS_THICKNESS
135           CALL TIMEAVE_CUMULATE( layers_HU_T, layers_HU, Nlayers,           CALL TIMEAVE_CUMULATE( layers_Hw_T, layers_HU, Nlayers,
136         &                          deltaTclock, bi, bj, myThid )
137             CALL TIMEAVE_CUMULATE( layers_PIw_T, layers_PIw, Nlayers,
138         &                          deltaTclock, bi, bj, myThid )
139             CALL TIMEAVE_CUMULATE( layers_U_T, layers_U, Nlayers,
140       &                          deltaTclock, bi, bj, myThid )       &                          deltaTclock, bi, bj, myThid )
141  #endif /* LAYERS_THICKNESS */  #endif /* LAYERS_THICKNESS */
142  #endif /* LAYERS_UFLUX */  #endif /* LAYERS_UFLUX */
143  #ifdef LAYERS_VFLUX  #ifdef LAYERS_VFLUX
144           CALL TIMEAVE_CUMULATE( layers_VFlux_T, layers_VFlux, Nlayers,           CALL TIMEAVE_CUMULATE( layers_VH_T, layers_VFlux, Nlayers,
145       &                          deltaTclock, bi, bj, myThid )       &                          deltaTclock, bi, bj, myThid )
146  #ifdef LAYERS_THICKNESS  #ifdef LAYERS_THICKNESS
147           CALL TIMEAVE_CUMULATE( layers_HV_T, layers_HV, Nlayers,           CALL TIMEAVE_CUMULATE( layers_Hs_T, layers_HV, Nlayers,
148         &                          deltaTclock, bi, bj, myThid )
149             CALL TIMEAVE_CUMULATE( layers_PIs_T, layers_PIs, Nlayers,
150         &                          deltaTclock, bi, bj, myThid )
151             CALL TIMEAVE_CUMULATE( layers_V_T, layers_V, Nlayers,
152       &                          deltaTclock, bi, bj, myThid )       &                          deltaTclock, bi, bj, myThid )
153  #endif /* LAYERS_THICKNESS */  #endif /* LAYERS_THICKNESS */
154  #endif /* LAYERS_VFLUX */  #endif /* LAYERS_VFLUX */
155    
156           layers_TimeAve(bi,bj)=layers_TimeAve(bi,bj)+deltaTclock  #ifdef LAYERS_PRHO_REF
157          IF (layers_num(iLa) .EQ. 3) THEN
158             CALL TIMEAVE_CUMULATE( prho_tave, prho, Nr,
159         &                          deltaTclock, bi, bj, myThid )
160        ENDIF        ENDIF
161  #endif /* ALLOW_TIMEAVE */  #endif /* LAYERS_PRHO_REF */
162    
163             layers_TimeAve(bi,bj)=layers_TimeAve(bi,bj)+deltaTclock
164    
165  C --- End bi,bj loop  C --- End bi,bj loop
166        ENDDO        ENDDO
167        ENDDO        ENDDO
168          ENDIF !IF ( layers_taveFreq.GT.0. ) THEN
169          ENDIF !IF ( iLa.EQ.1 ) THEN
170    #endif /* ALLOW_TIMEAVE */
171    
172    #endif /* ALLOW_LAYERS_OUTPUT */
173    
174          ENDDO !DO iLa=1,layers_maxNum
175          
176  #endif /* ALLOW_LAYERS */  #endif /* ALLOW_LAYERS */
177    
178        RETURN        RETURN
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22