/[MITgcm]/MITgcm/compare01/src/initialise.F
ViewVC logotype

Contents of /MITgcm/compare01/src/initialise.F

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

Revision 1.4 - (show annotations) (download)
Fri Feb 2 21:04:46 2001 UTC (23 years, 3 months ago) by adcroft
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +1 -1 lines
FILE REMOVED 
Merged changes from branch "branch-atmos-merge" into MAIN (checkpoint34)
 - substantial modifications to algorithm sequence (dynamics.F)
 - packaged OBCS, Shapiro filter, Zonal filter, Atmospheric Physics
1 C $Id: initialise.F,v 1.3 1998/06/12 19:33:31 cnh Exp $
2 #include "CPP_OPTIONS.h"
3 #include "CPP_MACROS.h"
4 #define SPHERICAL
5 #define _DB0 0.D0
6 C ================
7 C | initialise.F |
8 C ================
9 C Contents
10 C o INITIALISE Controlling procedure
11 C o INITALISE_CARTESIAN_GRID Initialise geometric terms in cartesian form
12 C o INITIALISE_CORIOLIS Initialise coriolis parameter
13 C o INITIALISE_DIAGS Initialise diagnostics.
14 C o INITIALISE_EOS Initialise equation of state operators
15 C o INITIALISE_FIELDS Set initial values of model fields
16 C o INITIALISE_GRIDDING Control setting up of model grid
17 C o INITIALISE_MASKS Set land-water masks
18 C o INITIALISE_METRIC Set curvilinear grid metric term operators
19 C o INITIALISE_OPERATORS Set finite difference operators
20 C o INITIALISE_S Set initial salinity field
21 C o INITIALISE_SPHERICAL_POLAR Initialise geometric terms in spherical polar form
22 C o INITIALISE_T Set inital temperature field
23 C o INITIALISE_U Set inital U velocity field
24 C o INITIALISE_V Set inital V velocity field
25
26 C================================================================================
27 C Procedure name: INITIALISE |
28 C Function: Controls initialisation of model fields and discrete |
29 C operators. |
30 C Comments: |
31 C================================================================================
32 CStartofinterface
33 SUBROUTINE INITIALISE (
34 O U, V, W, T, S, PH, PS )
35 IMPLICIT NONE
36 C /-------------------------------------------------------------------------\
37 C | Global variable declarations |
38 C \-------------------------------------------------------------------------/
39 #include "SIZE.h"
40 #include "PARAMS.h"
41 #include "STRINGS.h"
42 #include "OPERATORS.h"
43 #include "GRID.h"
44 #include "OLDG.h"
45 #include "MASKS.h"
46 #include "CG2DA.h"
47 #include "CG2DZ.h"
48 #include "AJAINF.h"
49 #include "FORCING.h"
50 C /-------------------------------------------------------------------------\
51 C | Routine argument declarations |
52 C |=========================================================================|
53 C | U, V, W - X,Y,Z Velocity ( m/s, m/s, Pa/s ). |
54 C | T - Potential temperature (oC). |
55 C | S - Salinity (ppt). |
56 C | PH - Hydrostatic pressure perturbation (m). |
57 C | PS - Surface pressure (m). |
58 C \-------------------------------------------------------------------------/
59 REAL U (_I3(nz,nx,ny))
60 REAL V (_I3(nz,nx,ny))
61 REAL W (_I3(nz,nx,ny))
62 REAL T (_I3(nz,nx,ny))
63 REAL S (_I3(nz,nx,ny))
64 REAL PH(_I3(nz,nx,ny))
65 REAL PS(nx, ny )
66 CEndofinteface
67 CALL INITIALISE_SPHERICAL_POLAR_GRID
68 CALL INITIALISE_CORIOLIS
69 CALL INITIALISE_FIELDS ( U, V, W, T, S, PH, PS )
70 C CALL INITIALISE_EOS
71 CALL INITIALISE_DIAGS
72 C
73 RETURN
74 END
75
76 C================================================================================
77 C Procedure name: INITIALISE_SPHERICAL_POLAR_GRID |
78 C Function: Iniitialise volumes and areas for a spherical polar |
79 C latitude-longitude grid. Basin geometry is specified in |
80 C terms of a southern boundary latitude (at V(j=1) ) and |
81 C latitudinal and longitudinal grid spacing in degrees. |
82 C Comments: |
83 C================================================================================
84 SUBROUTINE INITIALISE_SPHERICAL_POLAR_GRID
85 IMPLICIT NONE
86 C /-------------------------------------------------------------------------\
87 C | Global declarations |
88 C \-------------------------------------------------------------------------/
89 #include "SIZE.h"
90 #include "PARAMS.h"
91 #include "STRINGS.h"
92 #include "OPERATORS.h"
93 #include "GRID.h"
94 #include "CG2DA.h"
95 #include "CG2DZ.h"
96 #include "MASKS.h"
97 #include "FORCING.h"
98 C /-------------------------------------------------------------------------\
99 C | Local declarations |
100 C |=========================================================================|
101 C | sbLat - Latitude (in degress ) of southern boundary |
102 C | -90.0 = South Pole |
103 C | 0.0 = Equator |
104 C | +90.0 = North Pole |
105 C | delX - X grid spacing in degress longitude |
106 C | delY - Y grid spacing in degress latitude |
107 C | latSouth- Holds active latitude (degrees) |
108 C | latNorth- Holds active latitude (degrees) |
109 C | longWest- Holds active longitude (degrees) |
110 C | longEast- Holds active longitude (degrees) |
111 C | I,J,K - Loop counters |
112 C \-------------------------------------------------------------------------/
113 REAL sbLat
114 REAL delX
115 REAL delY
116 REAL latSouth, latNorth
117 REAL longWest, longEast
118 INTEGER I, J, K
119 REAL dist
120 REAL tmpVal
121 REAL tmpMsk( _I3(nz,0:nx+1,0:ny+1) )
122 REAL aC( 1:nx,1:ny )
123 LOGICAL landHere, waterHere
124 LOGICAL landSouth, waterSouth
125 LOGICAL landWest, waterWest
126 LOGICAL waterAbove, landAbove
127 LOGICAL TOP_LAYER, BOTTOM_LAYER
128 C
129 delX = dTheta
130 delY = dPhi
131 sbLat = phiMin
132 C
133 C
134 rZero = radiusEarth
135 #ifdef SPHERICAL
136 rPVolHat = 0.
137 DO k=1, nz
138 latNorth = sbLat
139 DO j =1, ny
140 latSouth = latNorth
141 latNorth = latNorth+delY
142 longEast = 0.
143 DO i =1, nx
144 longWest = longEast
145 longEast = longWest+delX
146 XA( _I3(k,i,j) ) = rZero*delY*deg2rad*delps(k)
147 YA( _I3(k,i,j) ) = rZero*delX*deg2rad*COS(latSouth*deg2rad)*delps(k)
148 ZA( _I3(k,i,j) ) = rZero*delX*deg2rad
149 & *(SIN(latNorth*deg2rad)-SIN(latSouth*deg2rad))
150 & *rZero
151 rPVol( _I3(k,i,j) )= 1.D0/( ZA( _I3(k,i,j) )*delps(k) )
152 rPVolHat(i,j) = rPVolHat(i,j) + ZA( _I3(k,i,j) )*delps(k)
153 dist = rZero*delX*deg2rad*COS( 0.5D0*deg2rad*(latNorth+latSouth) )
154 uDdxOp( _I3(k,i,j) ) = 1.D0/( dist )
155 pDdxOp( _I3(k,i,j) ) = 1.D0/( dist )
156 dist = rZero*delX*deg2rad*COS( deg2rad*latSouth )
157 vDdxOp( _I3(k,i,j) ) = 1.D0/( dist )
158 dist = rZero*delY*deg2rad
159 uDdyOp( _I3(k,i,j) ) = 1.D0/( dist )
160 vDdyOp( _I3(k,i,j) ) = 1.D0/( dist )
161 pDdyOp( _I3(k,i,j) ) = 1.D0/( dist )
162 C fCorV( _I3(k,i,j) ) = 2.D0*omega*SIN(latSouth*deg2rad)
163 fCorU( _I3(k,i,j) ) = 2.D0*omega*SIN( 0.5D0*deg2rad*(latNorth+latSouth) )
164 uTPhiOp( _I3(k,i,j) ) = TAN(0.5D0*deg2rad*(latNorth+latSouth) )
165 vTPhiOp( _I3(k,i,j) ) = TAN(deg2rad*(latSouth) )
166 uTPhiOp( _I3(k,i,j) ) = 0.
167 vTPhiOp( _I3(k,i,j) ) = 0.
168 ENDDO
169 ENDDO
170 ENDDO
171 fCorV = fCorV * _DB0
172 fCorV = 0.5D0*(fCorU+CSHIFT(fCorU,-1,2) )
173 C fCorU = 0.5*(fCorV+CSHIFT(fCorV,+1,2) )
174 fCorW = fCorU
175
176 C fCorU = 1.d-4
177 C fCorV = 1.d-4
178
179 WRITE(0,*) ' fcoru '
180 CALL PLOT_FIELD(fCorU(_I3(1,1:nx,1:ny)), Nx, Ny)
181
182 DO j=1, ny
183 DO i=1, nx
184 rPVolHat(i,j) = 1.D0/rPVolHat(i,j)
185 ENDDO
186 ENDDO
187 #endif
188 rDxatU = pDdxOp(_I3(1,:,:))
189 rDyatV = pDdyOp(_I3(1,:,:))
190 C
191 DO k=1, nz
192 DO j=1, ny
193 rUVol( _I3(k,1,j) ) = ZA( _I3(k,1 ,j) )*delps(k)
194 & +ZA( _I3(k,nx ,j) )*delps(k)
195 rUVol( _I3(k,1,j) ) = 2.D0/rUVol( _I3(k,1,j) )
196 DO i=2, nx
197 rUVol( _I3(k,i,j) ) = ZA( _I3(k,i ,j) )*delps(k)
198 & +ZA( _I3(k,i-1,j) )*delps(k)
199 rUVol( _I3(k,i,j) ) = 2.D0/rUVol( _I3(k,i,j) )
200 ENDDO
201 ENDDO
202 ENDDO
203 C
204 DO k=1, nz
205 DO j=2, ny
206 DO i=1, nx
207 rVVol( _I3(k,i,j) ) = ZA( _I3(k,i,j ) )*delps(k)
208 & +ZA( _I3(k,i,j-1) )*delps(k)
209 rVVol( _I3(k,i,j) ) = 2.D0/rVVol( _I3(k,i,j) )
210 ENDDO
211 ENDDO
212 ENDDO
213 DO k=1, nz
214 DO i=1, nx
215 rVVol( _I3(k,i,1) ) = ZA( _I3(k,i,1 ) )*delps(k)
216 & +ZA( _I3(k,i,ny ) )*delps(k)
217 rVVol( _I3(k,i,1) ) = 2.D0/rVVol( _I3(k,i,1) )
218 ENDDO
219 ENDDO
220 C
221 UMASK = WATER
222 VMASK = WATER
223 tmpMsk( _I3(:,1:Nx,1:Ny) ) = PMASK
224 tmpMsk( _I3(:,0 ,1:Ny) ) = tmpMsk( _I3(:,Nx,1:Ny) )
225 tmpMsk( _I3(:,Nx+1,1:Ny) ) = tmpMsk( _I3(:, 1,1:Ny) )
226 tmpMsk( _I3(:, :,0 ) ) = tmpMsk( _I3(:, :,Ny ) )
227 tmpMsk( _I3(:, :,Ny+1) ) = tmpMsk( _I3(:, :,1 ) )
228 DO k=1,nz
229 DO j=1,ny
230 DO i=1,nx
231 C du/dy boundary conditions.
232 WaterHere = tmpMsk( _I3(k,i ,j) )*tmpMsk( _I3(k,i-1,j) ) .GT. 0.5
233 LandHere = .NOT. WaterHere
234 WaterSouth = tmpMsk( _I3(k,i ,j-1) )*tmpMsk( _I3(k,i-1,j-1) ) .GT. 0.5
235 LandSouth = .NOT. WaterSouth
236 IF ( LandHere ) UMASK( _I3(k,i,j) ) = LAND
237 IF ( WaterHere .AND. LandSouth ) THEN
238 C uDdYop( _I3(k,i,j) ) = 0. ! Free slip southern boundaries
239 uDdYop( _I3(k,i,j) ) = +2.*uDdYop( _I3(k,i,j) ) ! No slip southern boundaries
240 ELSEIF ( LandHere .AND. WaterSouth ) THEN
241 C uDdYop( _I3(k,i,j) ) = 0. ! Free slip northern boundaries
242 uDdYop( _I3(k,i,j) ) = +2.*uDdYop( _I3(k,i,j) ) ! No slip northern boundaries
243 ELSEIF ( LandHere .AND. LandSouth ) THEN
244 uDdYop( _I3(k,i,j) ) = 0.
245 ENDIF
246 C dv/dx boundary conditions.
247 WaterHere = tmpMsk( _I3(k,i,j ) )*tmpMsk( _I3(k,i,j-1) ) .GT. 0.5
248 LandHere = .NOT. WaterHere
249 WaterWest = tmpMsk( _I3(k,i-1,j ) )*tmpMsk( _I3(k,i-1,j-1) ) .GT. 0.5
250 LandWest = .NOT. WaterWest
251 IF ( LandHere ) VMASK( _I3(k,i,j) ) = LAND
252 IF ( WaterHere .AND. LandWest ) THEN
253 C vDdXop( _I3(k,i,j) ) = 0. ! Free slip western boundaries
254 vDdXop( _I3(k,i,j) ) = +2.*vDdXop( _I3(k,i,j) ) ! No slip western boundaries
255 ELSEIF ( LandHere .AND. WaterWest ) THEN
256 C vDdXop( _I3(k,i,j) ) = 0. ! Free slip eastern boundaries
257 vDdXop( _I3(k,i,j) ) = +2.*vDdXop( _I3(k,i,j) ) ! No slip eastern boundaries
258 ELSEIF ( LandHere .AND. LandWest ) THEN
259 vDdXop( _I3(k,i,j) ) = 0.
260 ENDIF
261 C dp/dx and dp/dy boundary conditions
262 WaterHere = tmpMsk( _I3(k,i ,j ) ) .GT. 0.5
263 LandHere = .NOT. WaterHere
264 WaterSouth = tmpMsk( _I3(k,i ,j-1) ) .GT. 0.5
265 LandSouth = .NOT. WaterSouth
266 WaterWest = tmpMsk( _I3(k,i-1,j ) ) .GT. 0.5
267 LandWest = .NOT. WaterWest
268 IF ( LandHere ) THEN
269 pDdxOp( _I3(k,i,j) ) = 0.
270 pDdyOp( _I3(k,i,j) ) = 0.
271 ELSE
272 IF ( LandWest ) THEN
273 pDdxOp( _I3(k,i,j) ) = 0.
274 ENDIF
275 IF ( LandSouth ) THEN
276 pDdyOp( _I3(k,i,j) ) = 0.
277 ENDIF
278 ENDIF
279 ENDDO
280 ENDDO
281 C du/dz, dv/dz and dp/dz boundary conditions.
282 TOP_LAYER = K .EQ. 1
283 BOTTOM_LAYER = K .EQ. Nz
284 C Can skip first and last layers. These are handled in-line
285 C by special case code. The in-line TOP_LAYER, BOTTOM_LAYER
286 C code needs to be altered if the boundary conditions are altered.
287 C Only the land below case is coded up. This assumes that there
288 C are no "overhangs" in the topography.
289 IF ( .NOT. TOP_LAYER ) THEN
290 DO j=1,ny
291 DO i=1,nx
292 C du/dz = 0.
293 LandHere = tmpMsk( _I3(k ,i,j) )*tmpMsk( _I3(k ,i-1,j) ) .LT. 0.5
294 WaterAbove = tmpMsk( _I3(k-1,i,j) )*tmpMsk( _I3(k-1,i-1,j) ) .GT. 0.5
295 IF ( LandHere .AND. WaterAbove ) THEN
296 uDdzOp( _I3(k,i,j) ) = 0. ! Free Slip
297 uDdzOp( _I3(k,i,j) ) = rDzAtW(K)*2.D0 ! No Slip
298 ENDIF
299 C dv/dz = 0.
300 LandHere = tmpMsk( _I3(k ,i,j) )*tmpMsk( _I3(k ,i,j-1) ) .LT. 0.5
301 WaterAbove = tmpMsk( _I3(k-1,i,j) )*tmpMsk( _I3(k-1,i,j-1) ) .GT. 0.5
302 IF ( LandHere .AND. WaterAbove ) THEN
303 vDdzOp( _I3(k,i,j) ) = 0. ! Free Slip
304 vDdzOp( _I3(k,i,j) ) = rDzAtW(K)*2.D0 ! No Slip
305 ENDIF
306 C dt/dz = ds/dz = 0.
307 LandHere = tmpMsk( _I3(k ,i,j) ) .LT. 0.5
308 WaterHere = tmpMsk( _I3(k ,i,j) ) .GT. 0.5
309 WaterAbove = tmpMsk( _I3(k-1,i,j) ) .GT. 0.5
310 LandAbove = tmpMsk( _I3(k-1,i,j) ) .LT. 0.5
311 IF ( LandHere .AND. WaterAbove ) THEN
312 pDdzOp( _I3(k,i,j) ) = 0. ! Insulating
313 ENDIF
314 C Error overhanging shelf.
315 IF ( WaterHere .AND. LandAbove ) THEN
316 WRITE(scrUnitError,*) 'Land above water at I,J,K = ',I,J,K
317 initialisationError = .TRUE.
318 ENDIF
319 ENDDO
320 ENDDO
321 DO j=1,ny
322 DO i=1,nx
323 wMask(_I3(K,i,j)) = pMask(_I3(K,i,j))*pMask(_I3(K-1,i,j))
324 ENDDO
325 ENDDO
326 ENDIF
327 ENDDO
328 wMask(_I3(1,:,:)) = 0.
329 C
330 DO k=1, nz
331 DO j=1, ny
332 DO i=1, nx
333 rUVol( _I3(k,i,j) ) = rUVol( _I3(k,i,j) )*UMASK( _I3(k,i,j) )
334 rVVol( _I3(k,i,j) ) = rVVol( _I3(k,i,j) )*VMASK( _I3(k,i,j) )
335 rPVol( _I3(k,i,j) ) = rPVol( _I3(k,i,j) )*PMASK( _I3(k,i,j) )
336 ENDDO
337 ENDDO
338 ENDDO
339 XA = XA*UMASK
340 YA = YA*VMASK
341 ZA = ZA*PMASK
342 rPVolHat = rPVolHat*PMASK(_I3(1,:,:))
343 C
344 C Elliptic equation solver variables.
345 pC = 1.
346 pX = 0.
347 pY = 0.
348 aX = 0.
349 aY = 0.
350 DO K = 1, nz
351 aX(1:nx,1:ny) = aX(1:nx,1:ny)+XA(_I3(K,:,:))*pDdxOp(_I3(k,:,:))/RONIL
352 aY(1:nx,1:ny) = aY(1:nx,1:ny)+YA(_I3(K,:,:))*pDdyOp(_I3(k,:,:))/RONIL
353 ENDDO
354 aX(0 ,1:ny) = aX(nx,1:ny)
355 aX(nx+1,1:ny) = aX( 1,1:ny)
356 aY(1:nx,0 ) = aY(1:nx,ny+1)
357 aY(1:nx,ny+1) = aY(1:nx, 1)
358 aNorm2d = MAXVAL(ABS(aX))
359 tmpVal = MAXVAL(ABS(aY))
360 aNorm2d = MAX(aNorm2d,tmpVal)
361 DO J=1,ny
362 DO I=1,nx
363 pC(I,J) = -aX(I,J)-aX(I+1,J)
364 & -aY(I,J)-aY(I,J+1)
365 & -freeSurfFac*ZA(_I3(1,I,J))/DELT/DELT
366 aC(I,J) = pC(I,J)
367 IF ( pC(I,J) .NE. 0. ) pC(I,J) = 1.0/pC(I,J)
368 CcnhDebugStarts
369 C IF ( pC(I,J) .NE. 0. ) pC(I,J) = 1.D0
370 CcnhDebugEnds
371 ENDDO
372 ENDDO
373 aX = aX/aNorm2d
374 aY = aY/aNorm2d
375 aC = aC/aNorm2d
376 WHERE ( aC + CSHIFT(aC,-1,1) .NE. 0. )
377 pX(1:Nx,1:Ny) = -aX(1:Nx,1:Ny)/( (0.51*(aC + CSHIFT(aC,-1,1)))**2 )
378 ENDWHERE
379 WHERE ( aC + CSHIFT(aC,-1,2) .NE. 0. )
380 pY(1:Nx,1:Ny) = -aY(1:Nx,1:Ny)/( (0.51*(aC + CSHIFT(aC,-1,2)))**2 )
381 ENDWHERE
382 CcnhDebugStarts
383 pC = pC*aNorm2d
384 CcnhDebugEnds
385 C
386 CcnhDebugStarts
387 WRITE(0,*) ' aX = '
388 CALL PLOT_FIELD(aX(1:Nx,1:Ny), Nx, Ny)
389 WRITE(0,*) ' aY = '
390 CALL PLOT_FIELD(aY(1:Nx,1:Ny), Nx, Ny)
391 CcnhDebugEnds
392
393 RETURN
394 END
395 SUBROUTINE INITIALISE_CORIOLIS
396 IMPLICIT NONE
397 RETURN
398 END
399 C/-------------------------------------------------------------------\
400 C||| Procedure: INITIALISE_DIAGS |||
401 C|||===============================================================|||
402 C||| Function: Set up quantities that will be diagnosed. |||
403 C||| Comments: |||
404 C\-------------------------------------------------------------------/
405 CStartofinterface
406 SUBROUTINE INITIALISE_DIAGS
407 IMPLICIT NONE
408 C /--------------------------------------------------------------\
409 C | Global data |
410 C \--------------------------------------------------------------/
411 #include "SIZE.h"
412 #include "PARAMS.h"
413 #include "DIAGS.h"
414 CEndofinterface
415 C /--------------------------------------------------------------\
416 C | Routine arguments |
417 C |==============================================================|
418 C | ** NONE *** |
419 C \--------------------------------------------------------------/
420 C /--------------------------------------------------------------\
421 C | Local variables |
422 C |==============================================================|
423 C | I - Loop counter |
424 C \--------------------------------------------------------------/
425 INTEGER iErr
426 C /--------------------------------------------------------------\
427 C | Set basins. |
428 C \--------------------------------------------------------------/
429 C iErr = 0
430 C CALL DIAGS_ADD_BASIN('GLB',iErr)
431 C IF ( iErr .NE. 0 ) THEN
432 C WRITE(scrUnitError,*)
433 C & '*** Error: Basin GLB was not defined correctly.'
434 C initialisationError = .TRUE.
435 C ENDIF
436 C CALL DIAGS_ADD_BASIN('NA' ,iErr)
437 C IF ( iErr .NE. 0 ) THEN
438 C WRITE(scrUnitError,*)
439 C & '*** Error: Basin NA was not defined correctly.'
440 C initialisationError = .TRUE.
441 C ENDIF
442 C CALL DIAGS_ADD_BASIN('IO' ,iErr)
443 C IF ( iErr .NE. 0 ) THEN
444 C WRITE(scrUnitError,*)
445 C & '*** Error: Basin IO was not defined correctly.'
446 C initialisationError = .TRUE.
447 C ENDIF
448
449 CALL DIAGS_SHOW
450
451 RETURN
452 END
453
454 SUBROUTINE INITIALISE_FIELDS(
455 U U, V, W, T, S, PH, PS )
456 IMPLICIT NONE
457 C /-------------------------------------------------------------------------\
458 C | Global variable declarations |
459 C \-------------------------------------------------------------------------/
460 #include "SIZE.h"
461 #include "PARAMS.h"
462 #include "STRINGS.h"
463 #include "OPERATORS.h"
464 #include "GRID.h"
465 #include "OLDG.h"
466 #include "MASKS.h"
467 #include "CG2DA.h"
468 #include "CG2DZ.h"
469 #include "AJAINF.h"
470 #include "FORCING.h"
471 C /-------------------------------------------------------------------------\
472 C | Routine argument declarations |
473 C |=========================================================================|
474 C | U, V, W - X,Y,Z Velocity ( m/s, m/s, Pa/s ). |
475 C | T - Potential temperature (oC). |
476 C | S - Salinity (ppt). |
477 C | PH - Hydrostatic pressure perturbation (m). |
478 C | PS - Surface pressure (m). |
479 C \-------------------------------------------------------------------------/
480 REAL U (_I3(nz,nx,ny))
481 REAL V (_I3(nz,nx,ny))
482 REAL W (_I3(nz,nx,ny))
483 REAL T (_I3(nz,nx,ny))
484 REAL S (_I3(nz,nx,ny))
485 REAL PH(_I3(nz,nx,ny))
486 REAL PS(nx, ny )
487 CEndofinteface
488 C /-------------------------------------------------------------------------\
489 C | Local variables |
490 C |=========================================================================|
491 C | tmp Single precision I/O buffer |
492 C | iErr Error flag |
493 C | K Loop counter |
494 C \-------------------------------------------------------------------------/
495 Real*4 tmp(Nx,Ny,Nz)
496 INTEGER iErr
497 INTEGER K
498 IF ( PickupRun ) THEN
499 C /------------------------------------------------------------------------\
500 C | Run starts from previously saved model state. |
501 C |========================================================================|
502 C | This block handles "pickup" runs that are used to continue integrations|
503 C | over a series of computational jobs. |
504 C \------------------------------------------------------------------------/
505 iErr = 0
506 CcnhDEbugStarts
507 WRITE(0,*) ' Calling READ_DUMP'
508 CcnhDebugEnds
509 CALL READ_DUMP_3D(U ,Nx,Ny,Nz,'U.' ,puSuffix,iErr)
510 CcnhDebugStarts
511 WRITE(0,*) MAXVAL(U)
512 WRITE(0,*) MINVAL(U)
513 CcnhDEbugEnds
514 CALL READ_DUMP_3D(V ,Nx,Ny,Nz,'V.' ,puSuffix,iErr)
515 CALL READ_DUMP_3D(uAja ,Nx,Ny,Nz,'UAJA.' ,puSuffix,iErr)
516 CALL READ_DUMP_3D(vAja ,Nx,Ny,Nz,'VAJA.' ,puSuffix,iErr)
517 CALL READ_DUMP_3D(W ,Nx,Ny,Nz,'W.' ,puSuffix,iErr)
518 CcnhDebugStarts
519 CALL READ_DUMP_3D(T ,Nx,Ny,Nz,'T.' ,puSuffix,iErr)
520 CALL READ_DUMP_3D(S ,Nx,Ny,Nz,'S.' ,puSuffix,iErr)
521 CcnhDebugEnds
522 CALL READ_DUMP_2D(PS ,Nx,Ny, 'PS.' ,puSuffix,iErr)
523 CALL READ_DUMP_3D(PH ,Nx,Ny,Nz,'PH.' ,puSuffix,iErr)
524 CALL READ_DUMP_3D(UNM1 ,Nx,Ny,Nz,'UNM1.' ,puSuffix,iErr)
525 CALL READ_DUMP_3D(VNM1 ,Nx,Ny,Nz,'VNM1.' ,puSuffix,iErr)
526 CcnhDebugStarts
527 C Added cos. old code contains uBxBYNM1 etc which is not
528 C saved here.
529 UNM1 = 0.
530 VNM1 = 0.
531 CcnhDebugEnds
532 CALL READ_DUMP_2D(PSNM1,Nx,Ny, 'PSNM1.',puSuffix,iErr)
533 CALL READ_DUMP_3D(PHNM1,Nx,Ny,Nz,'PHNM1.',puSuffix,iErr)
534 CALL READ_DUMP_3D(GUNM1,Nx,Ny,Nz,'GUNM1.',puSuffix,iErr)
535 CALL READ_DUMP_3D(GVNM1,Nx,Ny,Nz,'GVNM1.',puSuffix,iErr)
536 CcnhDebugStarts
537 CALL READ_DUMP_3D(GTNM1,Nx,Ny,Nz,'GTNM1.',puSuffix,iErr)
538 CALL READ_DUMP_3D(GSNM1,Nx,Ny,Nz,'GSNM1.',puSuffix,iErr)
539 CcnhDebugEnds
540 IF ( iErr .NE. 0 ) THEN
541 C /----------------------------------------------------------------------\
542 C | Stop if there was error loading pickup data. |
543 C \----------------------------------------------------------------------/
544 CALL PRINT_ERRORS
545 initialisationError = .TRUE.
546 ENDIF
547 ELSE
548 C /------------------------------------------------------------------------\
549 C | Run starts from some externally generated state. |
550 C |========================================================================|
551 C | This block handles "startup" runs that typically start from t=0 and |
552 C | initialise the model fields from a climatology or with some predefined |
553 C | stratification and/or velcoity profile. |
554 C \------------------------------------------------------------------------/
555 C /------------------------------------------------------------------------\
556 C | Initialise velocity field. |
557 C |========================================================================|
558 C | Initial conditions set velocity field to zero. |
559 C | If non-zero velocity profile is used to initialise the model it needs |
560 C | to satisfy continuity to machine accuracy. |
561 C \------------------------------------------------------------------------/
562 U = 0.
563 V = 0.
564 W = 0.
565 C /------------------------------------------------------------------------\
566 C | Load potential temperaturea and salinity. |
567 C |========================================================================|
568 C | This example reads data derived from Levitus climatological average. |
569 C | Temperature is in degrees centigrade referenced to the surface. |
570 C | Salinity is in ppt. |
571 C \------------------------------------------------------------------------/
572 C /------------------------------------------------------------------------\
573 C | Temperature |
574 C \------------------------------------------------------------------------/
575 OPEN ( 11, FILE='LevCli.ptmp.sun.b', STATUS='old' , FORM='unformatted')
576 READ (11) tmp
577 CLOSE(11)
578 do K = 1, Nz
579 T( _I3(K,:,:) ) = tmp(:,:,K)
580 enddo
581 T = T*pMask
582 WRITE(0,*) 'T LOADED FROM FILE'
583 C /-----------------------------------------------------------------------\
584 C | Print character contour map of surface and deep temperature field. |
585 C \-----------------------------------------------------------------------/
586 CALL PLOT_FIELD(T(_I3(1,:,: )), Nx, Ny)
587 CALL PLOT_FIELD(T(_I3(20,:,:)), Nx, Ny)
588 C STOP
589 C /-----------------------------------------------------------------------\
590 C | Salinity |
591 C \-----------------------------------------------------------------------/
592 OPEN ( 11, FILE='LevCli.salt.sun.b', STATUS='old' , FORM='unformatted')
593 READ (11) tmp
594 CLOSE(11)
595 do K = 1, Nz
596 S( _I3(K,:,:) ) = tmp(:,:,K)
597 enddo
598 S = S*pMask
599 C /-----------------------------------------------------------------------\
600 C | Print character contour map of surface and deep salinity field. |
601 C \-----------------------------------------------------------------------/
602 C CALL PLOT_FIELD(S(_I3(1,:,: )), Nx, Ny)
603 C CALL PLOT_FIELD(S(_I3(20,:,:)), Nx, Ny)
604 ENDIF
605 C
606 RETURN
607 END
608 SUBROUTINE INITIALISE_EOS
609 IMPLICIT NONE
610 #include "SIZE.h"
611 #include "PARAMS.h"
612 #include "POLYEOS.h"
613 C sD - Salinity anomaly from reference
614 C tD - Temperature anomaly from reference
615 REAL sD
616 REAL tD
617 C I,J,K - Loop counters
618 INTEGER I, J, K
619 C Set polynomial parameters
620 OPEN( dUnit, FILE='polyeos_coeffs', STATUS='OLD' )
621 C READ ( dUnit, * ) tRef
622 C READ ( dUnit, * ) sRef
623 C READ ( dUnit, * ) C
624 C CLOSE( dUnit )
625 DO K = 1, Nz+1
626 READ(dUnit,*) tRef(K)
627 ENDDO
628 DO K = 1, Nz+1
629 READ(dUnit,*) sRef(K)
630 ENDDO
631 DO K = 1, Nz+1
632 DO I = 1, 9
633 READ(dUnit,*) C(K,I)
634 ENDDO
635 ENDDO
636 DO K = 1, Nz+1
637 READ(dUnit,*) sig0(K)
638 ENDDO
639 CLOSE(dUnit)
640 C Evaluate perturbation for T and S background state.
641 DO 10 K = 1, Nz
642 tD = ths(K)-tRef(K)
643 sD = (ssppt(K)-35.)/1000.-sRef(K)
644 rprmInit(K) = 1000.*(
645 & C(K,1)*tD +C(K,2)*sD
646 & +C(K,3)*tD**2+C(K,4)*tD*sD +C(K,5)*sD**2
647 & +C(K,6)*tD**3+C(K,7)*tD*sD**2+C(K,8)*sD*tD**2+C(K,9)*sD**3 )
648 10 CONTINUE
649
650 RETURN
651 END
  ViewVC Help
Powered by ViewVC 1.1.22