/[MITgcm]/MITgcm/pkg/seaice/seaice_init_varia.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_init_varia.F

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

revision 1.84 by mlosch, Mon Apr 28 15:17:32 2014 UTC revision 1.85 by heimbach, Mon May 12 22:08:54 2014 UTC
# Line 59  C     i,j,k,bi,bj :: Loop counters Line 59  C     i,j,k,bi,bj :: Loop counters
59         kSurface = 1         kSurface = 1
60        ENDIF        ENDIF
61    
62    C--   Initialize grid info
63          DO bj=myByLo(myThid),myByHi(myThid)
64           DO bi=myBxLo(myThid),myBxHi(myThid)
65            DO j=1-OLy,sNy+OLy
66             DO i=1-OLx,sNx+OLx
67              HEFFM(i,j,bi,bj)       = 0. _d 0
68             ENDDO
69            ENDDO
70            DO j=1-OLy,sNy+OLy
71             DO i=1-OLx,sNx+OLx
72              HEFFM(i,j,bi,bj)= 1. _d 0
73              IF (_hFacC(i,j,kSurface,bi,bj).eq.0.)
74         &         HEFFM(i,j,bi,bj)= 0. _d 0
75             ENDDO
76            ENDDO
77    #ifndef SEAICE_CGRID
78            DO j=1-OLy+1,sNy+OLy
79             DO i=1-OLx+1,sNx+OLx
80              UVM(i,j,bi,bj)=0. _d 0
81              mask_uice=HEFFM(i,j,  bi,bj)+HEFFM(i-1,j-1,bi,bj)
82         &             +HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j,  bi,bj)
83              IF(mask_uice.GT.3.5 _d 0) UVM(i,j,bi,bj)=1. _d 0
84             ENDDO
85            ENDDO
86    #endif /* SEAICE_CGRID */
87           ENDDO
88          ENDDO
89    
90    C     coefficients for metric terms
91          DO bj=myByLo(myThid),myByHi(myThid)
92           DO bi=myBxLo(myThid),myBxHi(myThid)
93    #ifdef SEAICE_CGRID
94            DO j=1-OLy,sNy+OLy
95             DO i=1-OLx,sNx+OLx
96              k1AtC(I,J,bi,bj) = 0.0 _d 0
97              k1AtZ(I,J,bi,bj) = 0.0 _d 0
98              k2AtC(I,J,bi,bj) = 0.0 _d 0
99              k2AtZ(I,J,bi,bj) = 0.0 _d 0
100             ENDDO
101            ENDDO
102            IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
103    C     This is the only case where tan(phi) is not zero. In this case
104    C     C and U points, and Z and V points have the same phi, so that we
105    C     only need a copy here. Do not use tan(YC) and tan(YG), because
106    C     these
107    C     can be the geographical coordinates and not the correct grid
108    C     coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
109             DO j=1-OLy,sNy+OLy
110              DO i=1-OLx,sNx+OLx
111               k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
112               k2AtZ(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
113              ENDDO
114             ENDDO
115            ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
116    C     compute metric term coefficients from finite difference
117    C     approximation
118             DO j=1-OLy,sNy+OLy
119              DO i=1-OLx,sNx+OLx-1
120               k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
121         &          * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
122         &          * _recip_dxF(I,J,bi,bj)
123              ENDDO
124             ENDDO
125             DO j=1-OLy,sNy+OLy
126              DO i=1-OLx+1,sNx+OLx
127               k1AtZ(I,J,bi,bj) = _recip_dyU(I,J,bi,bj)
128         &          * ( _dyC(I,J,bi,bj) - _dyC(I-1,J,bi,bj) )
129         &          * _recip_dxV(I,J,bi,bj)
130              ENDDO
131             ENDDO
132             DO j=1-OLy,sNy+OLy-1
133              DO i=1-OLx,sNx+OLx
134               k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
135         &          * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
136         &          * _recip_dyF(I,J,bi,bj)
137              ENDDO
138             ENDDO
139             DO j=1-OLy+1,sNy+OLy
140              DO i=1-OLx,sNx+OLx
141               k2AtZ(I,J,bi,bj) = _recip_dxV(I,J,bi,bj)
142         &          * ( _dxC(I,J,bi,bj) - _dxC(I,J-1,bi,bj) )
143         &          * _recip_dyU(I,J,bi,bj)
144              ENDDO
145             ENDDO
146            ENDIF
147    #else /* not SEAICE_CGRID */
148            DO j=1-OLy,sNy+OLy
149             DO i=1-OLx,sNx+OLx
150              k1AtC(I,J,bi,bj) = 0.0 _d 0
151              k1AtU(I,J,bi,bj) = 0.0 _d 0
152              k1AtV(I,J,bi,bj) = 0.0 _d 0
153              k2AtC(I,J,bi,bj) = 0.0 _d 0
154              k2AtU(I,J,bi,bj) = 0.0 _d 0
155              k2AtV(I,J,bi,bj) = 0.0 _d 0
156             ENDDO
157            ENDDO
158            IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
159    C     This is the only case where tan(phi) is not zero. In this case
160    C     C and U points, and Z and V points have the same phi, so that we
161    C     only need a copy here. Do not use tan(YC) and tan(YG), because
162    C     these
163    C     can be the geographical coordinates and not the correct grid
164    C     coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
165             DO j=1-OLy,sNy+OLy
166              DO i=1-OLx,sNx+OLx
167               k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
168               k2AtU(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
169               k2AtV(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
170              ENDDO
171             ENDDO
172            ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
173    C     compute metric term coefficients from finite difference
174    C     approximation
175             DO j=1-OLy,sNy+OLy
176              DO i=1-OLx,sNx+OLx-1
177               k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
178         &          * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
179         &          * _recip_dxF(I,J,bi,bj)
180              ENDDO
181             ENDDO
182             DO j=1-OLy,sNy+OLy
183              DO i=1-OLx+1,sNx+OLx
184               k1AtU(I,J,bi,bj) = _recip_dyG(I,J,bi,bj)
185         &          * ( _dyF(I,J,bi,bj) - _dyF(I-1,J,bi,bj) )
186         &          * _recip_dxC(I,J,bi,bj)
187              ENDDO
188             ENDDO
189             DO j=1-OLy,sNy+OLy
190              DO i=1-OLx,sNx+OLx-1
191               k1AtV(I,J,bi,bj) = _recip_dyC(I,J,bi,bj)
192         &          * ( _dyU(I+1,J,bi,bj) - _dyU(I,J,bi,bj) )
193         &          * _recip_dxG(I,J,bi,bj)
194              ENDDO
195             ENDDO
196             DO j=1-OLy,sNy+OLy-1
197              DO i=1-OLx,sNx+OLx
198               k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
199         &          * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
200         &          * _recip_dyF(I,J,bi,bj)
201              ENDDO
202             ENDDO
203             DO j=1-OLy,sNy+OLy-1
204              DO i=1-OLx,sNx+OLx
205               k2AtU(I,J,bi,bj) = _recip_dxC(I,J,bi,bj)
206         &          * ( _dxV(I,J+1,bi,bj) - _dxV(I,J,bi,bj) )
207         &          * _recip_dyG(I,J,bi,bj)
208              ENDDO
209             ENDDO
210             DO j=1-OLy+1,sNy+OLy
211              DO i=1-OLx,sNx+OLx
212               k2AtV(I,J,bi,bj) = _recip_dxG(I,J,bi,bj)
213         &          * ( _dxF(I,J,bi,bj) - _dxF(I,J-1,bi,bj) )
214         &          * _recip_dyC(I,J,bi,bj)
215              ENDDO
216             ENDDO
217            ENDIF
218    #endif /* not SEAICE_CGRID */
219           ENDDO
220          ENDDO
221    
222    #ifndef SEAICE_CGRID
223    C--   Choose a proxy level for geostrophic velocity,
224          DO bj=myByLo(myThid),myByHi(myThid)
225           DO bi=myBxLo(myThid),myBxHi(myThid)
226            DO j=1-OLy,sNy+OLy
227             DO i=1-OLx,sNx+OLx
228              KGEO(i,j,bi,bj)   = 0
229             ENDDO
230            ENDDO
231            DO j=1-OLy,sNy+OLy
232             DO i=1-OLx,sNx+OLx
233    #ifdef SEAICE_BICE_STRESS
234              KGEO(i,j,bi,bj) = 1
235    #else /* SEAICE_BICE_STRESS */
236              IF (klowc(i,j,bi,bj) .LT. 2) THEN
237               KGEO(i,j,bi,bj) = 1
238              ELSE
239               KGEO(i,j,bi,bj) = 2
240               DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND.
241         &          KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
242                  KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
243               ENDDO
244              ENDIF
245    #endif /* SEAICE_BICE_STRESS */
246             ENDDO
247            ENDDO
248           ENDDO
249          ENDDO
250    #endif /* SEAICE_CGRID */
251    
252    
253  C--   Initialise all variables in common blocks:  C--   Initialise all variables in common blocks:
254        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
255         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
Legend:
Removed from v.1.84  
changed lines
  Added in v.1.85
  ViewVC Help
Powered by ViewVC 1.1.22