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

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

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


Revision 1.3 - (hide annotations) (download)
Fri Jun 12 19:33:32 1998 UTC (26 years ago) by cnh
Branch: MAIN
Changes since 1.2: +124 -11 lines
Chages to make default setup correct for 4 degreee global comparison

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

  ViewVC Help
Powered by ViewVC 1.1.22