/[MITgcm]/MITgcm_contrib/rpa_layers/layers/layers_init_fixed.F
ViewVC logotype

Contents of /MITgcm_contrib/rpa_layers/layers/layers_init_fixed.F

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


Revision 1.1 - (show annotations) (download)
Tue Sep 15 19:16:53 2009 UTC (15 years, 10 months ago) by rpa
Branch: MAIN
importing layers package

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

  ViewVC Help
Powered by ViewVC 1.1.22