pkg/layers/layers_init_fixed.F

C $Header: /u/gcmpack/MITgcm/pkg/layers/layers_init_fixed.F,v 1.7 2012/09/19 18:48:18 gforget Exp $
C $Name:  $

#include "LAYERS_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      SUBROUTINE LAYERS_INIT_FIXED( myThid )

C ===================================================================
C     Initialize LAYERS variables that are kept fixed during the run.
C ===================================================================

      IMPLICIT NONE
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "LAYERS_SIZE.h"
#include "LAYERS.h"

C  INPUT/OUTPUT PARAMETERS:
C     myThid ::  my Thread Id number
      INTEGER myThid

C  LOCAL VARIABLES:
C     k         :: loop index
C     kk,kkinit :: fine grid loop indices
C     Zf        :: depth at cell boundaries
C     Zf        :: depth at cell centers
C     ZZf       :: depth at cell boundaries (fine grid)
C     ZZc       :: depth at cell centers (fine grid)
C     msgBuf    :: Informational/error message buffer
      INTEGER k,kk,kkinit,iLa
      _RL     Zf(Nr+1)
      _RL     Zc(Nr)
      _RL     ZZf(FineGridMax+1)
      _RL     ZZc(FineGridMax)

      CHARACTER*11   tmpName      
      CHARACTER*(MAX_LEN_MBUF) msgBuf

C     Functions ::
      INTEGER     ILNBLNK
      EXTERNAL    ILNBLNK
      
#ifdef ALLOW_MNC
#ifdef LAYERS_MNC
      IF (layers_MNC) THEN
        CALL LAYERS_MNC_INIT( myThid )
      ENDIF
#endif /* LAYERS_MNC */
#endif /* ALLOW_MNC */

C  Set up the vertical grid

C     for now, just use up the entire available array for ZZ
      NZZ = FineGridMax

C     Z and ZZ are INCREASING DOWNWARD!!!
C     Maybe this is dumb but it will work as long as we are consistent


C     Each dF cell is split into FineGridFact fine cells
C     Calculate dZZf on the fine grid
      kkinit = 1
      DO k=1,Nr
        DO kk=kkinit,kkinit+FineGridFact-1
          dZZf(kk) = dRf(k) / FineGridFact
        ENDDO
        kkinit = kkinit + FineGridFact
      ENDDO

C     find the depths
      Zf(1) = 0. _d 0
      Zc(1) = drC(1)
      DO k=2,Nr
        Zf(k) = Zf(k-1) + drF(k-1)
        Zc(k) = Zc(k-1) + drC(k)
      ENDDO
      Zf(Nr+1) = Zf(Nr) + drF(Nr)

C     create ZZ
      ZZf(1) = 0. _d 0
      ZZc(1) = 0.5 _d 0 * dZZf(1)

      DO kk=2,NZZ+1
            ZZf(kk) = ZZf(kk-1) + dZZf(kk-1)
            ZZc(kk-1) = 0.5 _d 0 * (ZZf(kk) + ZZf(kk-1))
      ENDDO

C     create the interpolating mapping arrays
      k = 1
      DO kk=1,NZZ
C       see if ZZc point is less than the top Zc point
        IF ( ZZc(kk) .LT. Zc(1) ) THEN
          MapIndex(kk) = 1
          MapFact(kk) = 1.0 _d 0
C       see if ZZc point is greater than the bottom Zc point
        ELSE IF ( (ZZc(kk) .GE. Zc(Nr)) .OR. (k .EQ. Nr) ) THEN
          MapIndex(kk) = Nr - 1
          MapFact(kk) = 0.0 _d 0
C       Otherwise we are somewhere in between and need to do interpolation)
        ELSE IF ( (ZZc(kk) .GE. Zc(k))
     &   .AND. (ZZc(kk) .LT. Zc(Nr)) ) THEN
C         Find the proper k value
          DO WHILE (ZZc(kk) .GE. Zc(k+1))
            k = k + 1
          ENDDO
C         If the loop stopped, that means Zc(k) <= ZZc(kk) < ZZc(k+1)
          MapIndex(kk) = k
          MapFact(kk) = 1.0 - (( ZZc(kk) - Zc(k) ) / drC(k+1))
        ELSE
C       This means there was a problem
          WRITE(msgBuf,'(A,I4,A,I4,A,1E14.6,A,2E14.6)')
     &     'S/R LAYERS_INIT_FIXED: kk=', kk, ', k=', k,
     &     ', ZZc(kk)=', ZZc(kk),' , Zc(k)=',Zc(k)
          CALL PRINT_ERROR( msgBuf, myThid )
          STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'
        ENDIF

C       See which grid box the point lies in
        IF ( ZZc(kk).LT.Zf(MapIndex(kk)+1) ) THEN
          CellIndex(kk) = MapIndex(kk)
        ELSE
          CellIndex(kk) = MapIndex(kk)+1
        ENDIF
      ENDDO

      IF ( debugLevel .GE. debLevB ) THEN
        CALL PRINT_MESSAGE( 'LAYERS_INIT_FIXED Debugging:',
     &             standardMessageUnit, SQUEEZE_RIGHT, myThid )
        DO kk=1,NZZ
          WRITE(msgBuf,'(A,1F6.1,A,I3,A,I3,A,I3,A,1F6.4,A,I3,A,I3)')
     &     '// ZZc=', ZZc(kk),
     &     ', MapIndex(',kk,')=',MapIndex(kk),
     &     ', MapFact(',kk,')=',MapFact(kk),
     &     ', CellIndex(',kk,')=',CellIndex(kk)
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )
        ENDDO
      ENDIF

C     Output the layer coordinates
      DO iLa=1,layers_maxNum
      IF ( layers_num(iLa).NE.0 ) THEN
      WRITE(tmpName,'(A7,I1,A3)') 'layers',iLa,layers_name(iLa)
      CALL WRITE_GLVEC_RL( tmpName, ' ', layers_G, 1+Nlayers,
     & -1, myThid )
      ENDIF
      ENDDO

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
        CALL LAYERS_DIAGNOSTICS_INIT( myThid )
      ENDIF
#endif

      RETURN
      END
1	gforget	1.8	C $Header: /u/gcmpack/MITgcm/pkg/layers/layers_init_fixed.F,v 1.7 2012/09/19 18:48:18 gforget Exp $
2	rpa	1.1	C $Name: $
3
4			#include "LAYERS_OPTIONS.h"
5
6			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
7
8			SUBROUTINE LAYERS_INIT_FIXED( myThid )
9
10			C ===================================================================
11			C Initialize LAYERS variables that are kept fixed during the run.
12			C ===================================================================
13
14			IMPLICIT NONE
15			#include "EEPARAMS.h"
16			#include "SIZE.h"
17			#include "PARAMS.h"
18			#include "GRID.h"
19			#include "LAYERS_SIZE.h"
20			#include "LAYERS.h"
21
22			C INPUT/OUTPUT PARAMETERS:
23			C myThid :: my Thread Id number
24			INTEGER myThid
25
26			C LOCAL VARIABLES:
27			C k :: loop index
28			C kk,kkinit :: fine grid loop indices
29			C Zf :: depth at cell boundaries
30			C Zf :: depth at cell centers
31			C ZZf :: depth at cell boundaries (fine grid)
32			C ZZc :: depth at cell centers (fine grid)
33	jmc	1.4	C msgBuf :: Informational/error message buffer
34	gforget	1.7	INTEGER k,kk,kkinit,iLa
35	rpa	1.1	_RL Zf(Nr+1)
36			_RL Zc(Nr)
37			_RL ZZf(FineGridMax+1)
38			_RL ZZc(FineGridMax)
39
40	gforget	1.7	CHARACTER*11 tmpName
41	rpa	1.1	CHARACTER*(MAX_LEN_MBUF) msgBuf
42
43	gforget	1.7	C Functions ::
44			INTEGER ILNBLNK
45			EXTERNAL ILNBLNK
46
47	rpa	1.1	#ifdef ALLOW_MNC
48			#ifdef LAYERS_MNC
49			IF (layers_MNC) THEN
50			CALL LAYERS_MNC_INIT( myThid )
51			ENDIF
52			#endif /* LAYERS_MNC */
53			#endif /* ALLOW_MNC */
54
55			C Set up the vertical grid
56
57			C for now, just use up the entire available array for ZZ
58			NZZ = FineGridMax
59
60			C Z and ZZ are INCREASING DOWNWARD!!!
61			C Maybe this is dumb but it will work as long as we are consistent
62
63
64			C Each dF cell is split into FineGridFact fine cells
65			C Calculate dZZf on the fine grid
66			kkinit = 1
67			DO k=1,Nr
68			DO kk=kkinit,kkinit+FineGridFact-1
69			dZZf(kk) = dRf(k) / FineGridFact
70			ENDDO
71			kkinit = kkinit + FineGridFact
72			ENDDO
73
74			C find the depths
75			Zf(1) = 0. _d 0
76			Zc(1) = drC(1)
77			DO k=2,Nr
78			Zf(k) = Zf(k-1) + drF(k-1)
79			Zc(k) = Zc(k-1) + drC(k)
80			ENDDO
81			Zf(Nr+1) = Zf(Nr) + drF(Nr)
82
83			C create ZZ
84			ZZf(1) = 0. _d 0
85	rpa	1.3	ZZc(1) = 0.5 _d 0 * dZZf(1)
86	rpa	1.1
87	rpa	1.3	DO kk=2,NZZ+1
88	rpa	1.1	ZZf(kk) = ZZf(kk-1) + dZZf(kk-1)
89	rpa	1.3	ZZc(kk-1) = 0.5 _d 0 * (ZZf(kk) + ZZf(kk-1))
90	rpa	1.1	ENDDO
91
92			C create the interpolating mapping arrays
93			k = 1
94			DO kk=1,NZZ
95			C see if ZZc point is less than the top Zc point
96			IF ( ZZc(kk) .LT. Zc(1) ) THEN
97			MapIndex(kk) = 1
98			MapFact(kk) = 1.0 _d 0
99			C see if ZZc point is greater than the bottom Zc point
100			ELSE IF ( (ZZc(kk) .GE. Zc(Nr)) .OR. (k .EQ. Nr) ) THEN
101			MapIndex(kk) = Nr - 1
102			MapFact(kk) = 0.0 _d 0
103			C Otherwise we are somewhere in between and need to do interpolation)
104			ELSE IF ( (ZZc(kk) .GE. Zc(k))
105			& .AND. (ZZc(kk) .LT. Zc(Nr)) ) THEN
106			C Find the proper k value
107			DO WHILE (ZZc(kk) .GE. Zc(k+1))
108			k = k + 1
109			ENDDO
110			C If the loop stopped, that means Zc(k) <= ZZc(kk) < ZZc(k+1)
111			MapIndex(kk) = k
112			MapFact(kk) = 1.0 - (( ZZc(kk) - Zc(k) ) / drC(k+1))
113			ELSE
114			C This means there was a problem
115			WRITE(msgBuf,'(A,I4,A,I4,A,1E14.6,A,2E14.6)')
116			& 'S/R LAYERS_INIT_FIXED: kk=', kk, ', k=', k,
117			& ', ZZc(kk)=', ZZc(kk),' , Zc(k)=',Zc(k)
118			CALL PRINT_ERROR( msgBuf, myThid )
119			STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'
120	jmc	1.4	ENDIF
121	rpa	1.1
122			C See which grid box the point lies in
123	jmc	1.2	IF ( ZZc(kk).LT.Zf(MapIndex(kk)+1) ) THEN
124	rpa	1.1	CellIndex(kk) = MapIndex(kk)
125			ELSE
126			CellIndex(kk) = MapIndex(kk)+1
127	jmc	1.4	ENDIF
128	rpa	1.1	ENDDO
129
130	jmc	1.5	IF ( debugLevel .GE. debLevB ) THEN
131	rpa	1.1	CALL PRINT_MESSAGE( 'LAYERS_INIT_FIXED Debugging:',
132	jmc	1.5	& standardMessageUnit, SQUEEZE_RIGHT, myThid )
133	rpa	1.1	DO kk=1,NZZ
134			WRITE(msgBuf,'(A,1F6.1,A,I3,A,I3,A,I3,A,1F6.4,A,I3,A,I3)')
135			& '// ZZc=', ZZc(kk),
136			& ', MapIndex(',kk,')=',MapIndex(kk),
137			& ', MapFact(',kk,')=',MapFact(kk),
138			& ', CellIndex(',kk,')=',CellIndex(kk)
139			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
140	jmc	1.5	& SQUEEZE_RIGHT, myThid )
141	rpa	1.1	ENDDO
142	jmc	1.4	ENDIF
143	rpa	1.1
144	rpa	1.3	C Output the layer coordinates
145	gforget	1.7	DO iLa=1,layers_maxNum
146			IF ( layers_num(iLa).NE.0 ) THEN
147			WRITE(tmpName,'(A7,I1,A3)') 'layers',iLa,layers_name(iLa)
148	rpa	1.9	CALL WRITE_GLVEC_RL( tmpName, ' ', layers_G, 1+Nlayers,
149	rpa	1.3	& -1, myThid )
150	gforget	1.7	ENDIF
151			ENDDO
152	rpa	1.3
153	gforget	1.6	#ifdef ALLOW_DIAGNOSTICS
154			IF ( useDiagnostics ) THEN
155			CALL LAYERS_DIAGNOSTICS_INIT( myThid )
156			ENDIF
157			#endif
158
159	rpa	1.1	RETURN
160			END