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

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

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


Revision 1.4 - (show annotations) (download)
Mon Jun 15 05:13:54 1998 UTC (26 years ago) by cnh
Branch: MAIN
CVS Tags: checkpoint7, checkpoint9, checkpoint8
Branch point for: checkpoint7-4degree-ref
Changes since 1.3: +31 -10 lines
Fairly coplete 4 degree global intercomparison
setup.
 Includes changes to make convective adjustment and hydrostatic
pressure correct as well as IO for climatological datasets

1 C $Id: set_defaults.F,v 1.3 1998/06/12 19:33:32 cnh Exp $
2 #include "CPP_OPTIONS.h"
3 #include "CPP_MACROS.h"
4 C================================================================================
5 C Procedure name: SET_DEFAULTS |
6 C Function: Assign default values to model parameters. |
7 C Comments: |
8 C================================================================================
9
10 CStartofinterface
11 SUBROUTINE SET_DEFAULTS (
12 O U, V, W, T, S, PH, PS )
13 IMPLICIT NONE
14 C /-------------------------------------------------------------------------\
15 C | Global variable declarations |
16 C \-------------------------------------------------------------------------/
17 #include "SIZE.h"
18 #include "PARAMS.h"
19 #include "STRINGS.h"
20 #include "OPERATORS.h"
21 #include "GRID.h"
22 #include "OLDG.h"
23 #include "MASKS.h"
24 #include "CG2DA.h"
25 #include "CG2DZ.h"
26 #include "AJAINF.h"
27 #include "EPARAM.h"
28 #include "FORCING.h"
29 C /-------------------------------------------------------------------------\
30 C | Routine argument declarations |
31 C |=========================================================================|
32 C | U, V, W - X,Y,Z Velocity ( m/s, m/s, Pa/s ). |
33 C | T - Potential temperature (oC). |
34 C | S - Salinity (ppt). |
35 C | PH - Hydrostatic pressure perturbation (m). |
36 C | PS - Surface pressure (m). |
37 C \-------------------------------------------------------------------------/
38 REAL U (_I3(nz,nx,ny))
39 REAL V (_I3(nz,nx,ny))
40 REAL W (_I3(nz,nx,ny))
41 REAL T (_I3(nz,nx,ny))
42 REAL S (_I3(nz,nx,ny))
43 REAL PH(_I3(nz,nx,ny))
44 REAL PS(nx, ny )
45
46 CEndofinteface
47
48 C /-------------------------------------------------------------------------\
49 C | Local variable declarations |
50 C |=========================================================================|
51 C | DELTAZ - Layer thickness (m). |
52 C | torry - Reference wind stress (N). |
53 C | y - Distance accumulator (m). |
54 C | ySize - Basin extent in y (m). |
55 C | pBeta - Planetary beta (s^-1.m^-1). |
56 C | I,J,K - Loop counters. |
57 C \-------------------------------------------------------------------------/
58 REAL DELTAZ(nz)
59 REAL torry
60 REAL y
61 REAL ySize
62 REAL pBeta
63 REAL Q
64 REAL DELTAF
65 REAL Fref
66 REAL fWall
67 REAL Beta
68 REAL r
69 REAL Kz
70 REAL Dep
71 REAL Rho0
72 REAL dist
73 REAL RMAX
74 REAL RSQ
75 REAL tau
76 REAL XCENTER
77 REAL YCENTER
78 REAL mx
79 REAL mn
80 REAL HEATFLAG
81 REAL HEATSUM
82 REAL RANDN
83 INTEGER I
84 INTEGER J
85 INTEGER K
86
87 REAL TAUMAX, DISTY, LY
88 Real*4 tmpXY(Nx,Ny)
89 Real*4 tmpXYZ( _I3(Nz,Nx,Ny) )
90 C
91 C /-------------------------------------------------------------------------\
92 C | Parameter settings for small double gyre beta-plane experiment. |
93 C \=========================================================================|
94 C
95 G = 9.81D0
96 RONIL = 999.8D0
97 CP = 3900.D0
98 CENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNORE
99
100 delT = 3600.D0
101 freeSurfFac = 1.D0
102 startTime = 0.
103 endTime = startTime+10.*delT
104 dPUFreq = 300.*delT
105
106
107 CSTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORE
108 phiMin = -80.
109 dPhi = 4.0
110 dTheta = 4.0
111 CENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNORE
112
113
114 a2MomXY = 5.0D5
115 a2MomZ = 1.D-3*G*RONIL*G*RONIL
116 C a2MomXY = 0.0D5
117 C a2MomZ = 0.D-3*G*RONIL*G*RONIL
118 a4MomXY = 0.D0
119 a2TempXY = 1.D3
120 a2TempZ = 3.D-5*G*RONIL*G*RONIL
121 a2SaltXY = 1.D3
122 a2SaltZ = 3.D-5*G*RONIL*G*RONIL
123 C a2SaltXY = 0.D3
124 C a2SaltZ = 0.D-5*G*RONIL*G*RONIL
125 C a2TempZ = 1.D-5*G*RONIL*G*RONIL
126 a4TempXy = 0.D0
127 a4SaltXy = 0.D0
128
129
130 CSTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORE
131 WRITE(20,*) ' delT',delT
132 WRITE(20,*) ' startTime',startTime
133 WRITE(20,*) ' endTime',endTime
134 WRITE(20,*) ' dPUFreq',dPUFreq
135 WRITE(20,*) ' a2MomXY',a2MomXY
136 WRITE(20,*) ' a2MomZ',a2MomZ
137 WRITE(20,*) ' a2TempXY',a2TempXY
138 WRITE(20,*) ' a2TempZ',a2TempZ
139 C********* GMGS *******************
140 C Ratio of time step
141 asyncfac = 30.0
142 C Scale dependent mixing scheme parameters
143 difref = a2tempXY
144 difmin = 0.1 * a2tempXY
145 difhmx = 2.5 * a2tempXY
146 coef = 1.8E9
147 dbotm = 1000.
148 C NOTE: MOD( INT(dslpcalcfreq), INT(dkcalcfreq) ) must = 0
149 dslpcalcfreq = 2. * delT
150 dkcalcfreq = 12. * delT
151 C Threashold for mix layer base
152 MixLayrDensGrad = 0.05
153 C Slope factors
154 slpmax1=5.e-3
155 slpmax2=1.e-2
156 slpmax3=1.e-3
157 small=1.e-30
158 dthin=50.
159 C Minimum density stratification
160 MinSigStrat= -1.e-8
161 c epsl = Av/Ah in Redis scheme
162 epsl = 0.
163 trDifScheme = gmconKh
164 verticalProfKh = .FALSE.
165
166 divHmax = 1.D+300
167 divHMin = 0.
168 rDelt = 1./DELT
169 numberOfTimeSteps = NINT((endTime-startTime)/delT)
170 WRITE(20,*) ' numberOfTimeSteps',numberOfTimeSteps
171 currentTime = startTime
172 nIter = NINT(startTime/delt)
173 WRITE(20,*) ' nIter',nIter
174 WRITE(puSuffix,'(I10.10)') nIter
175 PickupRun = .FALSE.
176 IF ( startTime .NE. 0. ) pickupRun = .TRUE.
177 C State variables
178 gsNM1 = 0.
179 gtNM1 = 0.
180 guNM1 = 0.
181 gvNM1 = 0.
182 U = 0.
183 V = 0.
184 W = 0.
185 S = 0.
186 PH = 0.
187 PS = 0.
188 PSNM1 = 0.
189 PHNM1 = 0.
190 uAja = 0.
191 vAja = 0.
192 CENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNORE
193
194
195 THS = 0.D0
196 THS(1) = 20.D0
197 THS(2) = 10.D0
198 THS(3) = 8.D0
199 THS(4) = 6.D0
200 SSPPT = 10.D0
201 T(_I3(1,:,:)) = 20.0D0
202 T(_I3(2,:,:)) = 10.0D0
203 T(_I3(3,:,:)) = 8.0D0
204 T(_I3(4,:,:)) = 6.0D0
205 C T(_I3(1,:,:)) = 1.D-12
206 C T(_I3(2,:,:)) = 0.
207 C T(_I3(3,:,:)) = 0.
208 C T(_I3(4,:,:)) = 0.
209 C T(_I3(1,:,:)) = 23.0D0
210 C T(_I3(2,:,:)) = 22.0D0
211 C T(_I3(3,:,:)) = 21.0D0
212 C T(_I3(4,:,:)) = 20.0D0
213 S = 10.0D0
214 THS( 1) = 16.0D0
215 THS( 2) = 15.2D0
216 THS( 3) = 14.5D0
217 THS( 4) = 13.9D0
218 THS( 5) = 13.3D0
219 THS( 6) = 12.4D0
220 THS( 7) = 11.3D0
221 THS( 8) = 9.9D0
222 THS( 9) = 8.4D0
223 THS(10) = 6.7D0
224 THS(11) = 5.2D0
225 THS(12) = 3.8D0
226 THS(13) = 2.9D0
227 THS(14) = 2.3D0
228 THS(15) = 1.8D0
229 THS(16) = 1.5D0
230 THS(17) = 1.1D0
231 THS(18) = 0.8D0
232 THS(19) = 0.66D0
233 THS(20) = 0.63D0
234 T(_I3( 1,:,:)) = THS( 1)
235 T(_I3( 2,:,:)) = THS( 2)
236 T(_I3( 3,:,:)) = THS( 3)
237 T(_I3( 4,:,:)) = THS( 4)
238 T(_I3( 5,:,:)) = THS( 5)
239 T(_I3( 6,:,:)) = THS( 6)
240 T(_I3( 7,:,:)) = THS( 7)
241 T(_I3( 8,:,:)) = THS( 8)
242 T(_I3( 9,:,:)) = THS( 9)
243 T(_I3(10,:,:)) = THS(10)
244 T(_I3(11,:,:)) = THS(11)
245 T(_I3(12,:,:)) = THS(12)
246 T(_I3(13,:,:)) = THS(13)
247 T(_I3(14,:,:)) = THS(14)
248 T(_I3(15,:,:)) = THS(15)
249 T(_I3(16,:,:)) = THS(16)
250 T(_I3(17,:,:)) = THS(17)
251 T(_I3(18,:,:)) = THS(18)
252 T(_I3(19,:,:)) = THS(19)
253 T(_I3(20,:,:)) = THS(20)
254
255 SSPPT( 1) = 34.65
256 SSPPT( 2) = 34.75
257 SSPPT( 3) = 34.82
258 SSPPT( 4) = 34.87
259 SSPPT( 5) = 34.90
260 SSPPT( 6) = 34.90
261 SSPPT( 7) = 34.86
262 SSPPT( 8) = 34.78
263 SSPPT( 9) = 34.69
264 SSPPT(10) = 34.60
265 SSPPT(11) = 34.58
266 SSPPT(12) = 34.62
267 SSPPT(13) = 34.68
268 SSPPT(14) = 34.72
269 SSPPT(15) = 34.73
270 SSPPT(16) = 34.74
271 SSPPT(17) = 34.73
272 SSPPT(18) = 34.73
273 SSPPT(19) = 34.72
274 SSPPT(20) = 34.72
275 S(_I3( 1,:,:)) = SSPPT( 1)
276 S(_I3( 2,:,:)) = SSPPT( 2)
277 S(_I3( 3,:,:)) = SSPPT( 3)
278 S(_I3( 4,:,:)) = SSPPT( 4)
279 S(_I3( 5,:,:)) = SSPPT( 5)
280 S(_I3( 6,:,:)) = SSPPT( 6)
281 S(_I3( 7,:,:)) = SSPPT( 7)
282 S(_I3( 8,:,:)) = SSPPT( 8)
283 S(_I3( 9,:,:)) = SSPPT( 9)
284 S(_I3(10,:,:)) = SSPPT(10)
285 S(_I3(11,:,:)) = SSPPT(11)
286 S(_I3(12,:,:)) = SSPPT(12)
287 S(_I3(13,:,:)) = SSPPT(13)
288 S(_I3(14,:,:)) = SSPPT(14)
289 S(_I3(15,:,:)) = SSPPT(15)
290 S(_I3(16,:,:)) = SSPPT(16)
291 S(_I3(17,:,:)) = SSPPT(17)
292 S(_I3(18,:,:)) = SSPPT(18)
293 S(_I3(19,:,:)) = SSPPT(19)
294 S(_I3(20,:,:)) = SSPPT(20)
295 DELTAZ( 1) = 50.D0
296 DELTAZ( 2) = 50.D0
297 DELTAZ( 3) = 55.D0
298 DELTAZ( 4) = 60.D0
299 DELTAZ( 5) = 65.D0
300 DELTAZ( 6) = 70.D0
301 DELTAZ( 7) = 80.D0
302 DELTAZ( 8) = 95.D0
303 DELTAZ( 9) = 120.D0
304 DELTAZ(10) = 155.D0
305 DELTAZ(11) = 200.D0
306 DELTAZ(12) = 260.D0
307 DELTAZ(13) = 320.D0
308 DELTAZ(14) = 400.D0
309 DELTAZ(15) = 480.D0
310 DELTAZ(16) = 570.D0
311 DELTAZ(17) = 655.D0
312 DELTAZ(18) = 725.D0
313 DELTAZ(19) = 775.D0
314 DELTAZ(20) = 815.D0
315
316 C Geometry
317 DM = 20. D3
318 DM = 4.
319 WRITE(20,*) ' DM',DM
320
321
322 CBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNORE
323 rPvolHat = 0.
324 DO K = 1, nz
325 DELPS(K) = G*RONIL*DELTAZ(K)
326 XA(_I3(K,:,:)) = DELPS(K)*DM
327 YA(_I3(K,:,:)) = DELPS(K)*DM
328 rUVol(_I3(K,:,:)) = 1.D0/DM/DM/DELPS(K)
329 rVVol(_I3(K,:,:)) = 1.D0/DM/DM/DELPS(K)
330 rPVol(_I3(K,:,:)) = 1.D0/DM/DM/DELPS(K)
331 rPvolHat = rPvolHat+DM*DM*DELPS(K)
332 IF ( K .EQ. 1 ) THEN
333 uDdzOp(_I3(K,:,:)) = 1.D0/DELPS(K)
334 vDdzOp(_I3(K,:,:)) = 1.D0/DELPS(K)
335 pDdzOp(_I3(K,:,:)) = 1.D0/DELPS(K)
336 ELSE
337 uDdzOp(_I3(K,:,:)) = 2.D0/(DELPS(K)+DELPS(K-1))
338 vDdzOp(_I3(K,:,:)) = 2.D0/(DELPS(K)+DELPS(K-1))
339 pDdzOp(_I3(K,:,:)) = 2.D0/(DELPS(K)+DELPS(K-1))
340 ENDIF
341 ! begin *** Redi Operators
342 rDZatP(K)=1.D0/DELPS(K)
343 if (K.eq.1) then
344 rDZatW(K)=1.D0/(DELPS(K))
345 else
346 rDZatW(K)=2.D0/(DELPS(K)+DELPS(K-1))
347 endif
348 ! end ***** Redi Operators
349 ENDDO
350 ! begin *** Redi Operators
351 rDXatU=1.d0/DM
352 rDYatV=1.d0/DM
353 ! end ***** Redi Operators
354
355 ZA = DM*DM
356 rPVolHat = 1.D0/rPvolHat
357 uDdxOp = 1.D0/DM
358 uDdyOp = 1.D0/DM
359 vDdxOp = uDdxOp
360 vDdyOp = uDdyOp
361 pDdxOp = uDdxOp
362 pDdyOp = uDdyOp
363 C Physical parameters
364 uTphiOp = 0.
365 vTphiOp = 0.
366 fCorW = 0.
367 pBeta = 1.D-11
368 CENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNORE
369
370
371 C radial BETA code begins - to use comment out 'NOBETA CODE' following
372 XCENTER = (Nx + 1)/2
373 WRITE(0,*) ' XCENTER: ', XCENTER
374 YCENTER = (Ny + 1)/2
375 WRITE(0,*) ' YCENTER: ', YCENTER
376 RMAX = 54.0
377 WRITE(20,*) ' RMAX: ',RMAX
378 WRITE(20,*) ' actual radius (m): ',(RMAX*DM)
379 fRef = 4.0
380 WRITE(20,*) ' fRef: ',fRef
381 fWall = 4.0
382 WRITE(20,*) ' fWall: ',fWall
383 Beta = (fWall - fRef)/(RMAX)
384 WRITE(20,*) ' Beta(s-1/gridpoint): ',Beta
385 WRITE(20,*) ' Beta(s-1/m): ',Beta/DM
386
387
388
389 C start NOBETACODE
390
391 fCorV(_I3(:,:,1)) = 1.D-4
392 Cd2
393 C pBeta = 0.
394 Cd2
395 write(20,*) 'NO RADIAL BETA'
396
397 fCorU(_I3(:,:,1)) = fCorV(_I3(:,:,1))+0.5*DM*pBeta
398 DO J = 2, ny
399 fCorV(_I3(:,:,J)) = fCorV(_I3(:,:,J-1))+DM*pBeta
400 fCorU(_I3(:,:,J)) = fCorU(_I3(:,:,J-1))+DM*pBeta
401 ENDDO
402
403 C end NOBETACODE
404
405
406 CBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNORE
407 sBeta = 0.
408 sBeta = 7.4E-4
409 WRITE(20,*) ' sBeta',sBeta
410 tAlpha = 2.D-4
411 C tAlpha = 0.D-4
412 WRITE(20,*) ' tAlpha',tAlpha
413
414
415 C Numerics parameters
416 abEps = 0.1
417 rAja = 9.887999800888535E-01
418 C rAja = 0.
419 dUnit = 10
420 scrUnit1 = 11
421 scrUnit2 = 12
422 scrUnitError = 13
423 InitialisationError = .FALSE.
424 IntegrationError = .FALSE.
425
426 C Topography
427 PMASK = WATER
428 C Enclose in a box
429 C Wall at X = Lx
430 PMASK(_I3(:,nx,:)) = LAND ! Land in east most cell
431 C Wall at Y = Ly
432 PMASK(_I3(:,:,ny)) = LAND ! Land in northern most cell
433 bathySource = binaryBathymetry
434 C Island
435 PMASK(_I3(Nz,1,24)) = LAND
436
437 C
438
439 MAX2DIT = 1000
440 TOLER2D = 1.D-13
441 freqCheckToler2d = 1
442
443 C-- Momentum equation forcing terms
444 tauMax = 0.1D0
445 lY = 0.
446 DO j=1,nY-1
447 lY = lY + DM
448 ENDDO
449 distY = 0.
450 DO j=1,nY
451 distY = distY+DM*0.5D0
452 DO i=1,nX
453 fv(i,j) = 0.D0
454 fu(i,j) = -tauMax*cos(2.D0*PI*distY/lY)
455 C fu(i,j) = -tauMax
456 fu(i,j) = tauMax*sin(PI*distY/lY)
457 fu(i,j) = fu(i,j)/(deltaZ(1)*ronil)
458 ENDDO
459 distY = distY+DM*0.5D0
460 ENDDO
461 fu(4,4) = fu(4,4)*0.917D0
462
463 OPEN ( 11, FILE='windx', STATUS='old' , FORM='unformatted')
464 READ (11) tmpXY
465 C Fu = tmpXY/DELPS(1)*G*0.1*uMask(_I3(1,:,:)) ! Dynes/cm**2 -> m/s**2
466 Fu = tmpXY/DELPS(1)*G*0.1
467 Fu = 0.1/DELPS(1)*G
468 CLOSE(11)
469
470 OPEN ( 11, FILE='windy', STATUS='old' , FORM='unformatted')
471 READ (11) tmpXY
472 C Fu = tmpXY/DELPS(1)*G*0.1*uMask(_I3(1,:,:)) ! Dynes/cm**2 -> m/s**2
473 Fv = tmpXY/DELPS(1)*G*0.
474 CLOSE(11)
475
476 C
477 CcnhDebugStarts
478 Cdbg WRITE(0,*) ' lY = ', lY
479 Cdbg WRITE(0,*) ' distY = ', distY
480 Cdbg WRITE(0,*) ' tauMax = ', tauMax
481 WRITE(0,*) ' fu '
482 CALL PLOT_FIELD(fu,Nx,Ny)
483 WRITE(0,*) ' fv '
484 CALL PLOT_FIELD(fv,Nx,Ny)
485 C STOP
486 CcnhDebugEnds
487
488 C T and S profiles
489 OPEN ( 11, FILE='LevCli.ptmp.sun.b', STATUS='old' , FORM='unformatted')
490 READ (11) tmpXYZ
491 Heat = tmpXYZ(_I3(1,:,:))
492 CLOSE(11)
493 C T and S profiles
494 OPEN ( 11, FILE='LevCli.salt.sun.b', STATUS='old' , FORM='unformatted')
495 READ (11) tmpXYZ
496 SSurf = tmpXYZ(_I3(1,:,:))
497 CLOSE(11)
498 C
499 END

  ViewVC Help
Powered by ViewVC 1.1.22