/[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.1 - (hide annotations) (download)
Mon May 25 20:21:05 1998 UTC (26 years, 1 month ago) by cnh
Branch: MAIN
CVS Tags: checkpoint4, checkpoint3, checkpoint5
Added version of compare01 reference code to repository.
Code committed is configured to produce same results as MITgcmUV

1 cnh 1.1 C $Id: set_defaults.F,v 1.16 1997/06/24 16:50:47 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     C
89     C /-------------------------------------------------------------------------\
90     C | Parameter settings for small double gyre beta-plane experiment. |
91     C \=========================================================================|
92     C
93     G = 9.81D0
94     RONIL = 999.8D0
95     CP = 3900.D0
96     CENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNORE
97    
98     delT = 1200.D0
99     freeSurfFac = 1.D0
100     startTime = 0.
101     endTime = startTime+10.*delT
102     dPUFreq = 300.*delT
103    
104    
105     CSTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORE
106     phiMin = 0.
107     dPhi = 4.0
108     dTheta = 4.0
109     CENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNORE
110    
111    
112     a2MomXY = 4.D2
113     a2MomZ = 1.D-2*G*RONIL*G*RONIL
114     a4MomXY = 0.D0
115     a2TempXY = 4.D2
116     a2TempZ = 1.D-2*G*RONIL*G*RONIL
117     C a2TempZ = 1.D-5*G*RONIL*G*RONIL
118     a4TempXy = 0.D0
119    
120    
121     CSTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORESTARTIGNORE
122     WRITE(20,*) ' delT',delT
123     WRITE(20,*) ' startTime',startTime
124     WRITE(20,*) ' endTime',endTime
125     WRITE(20,*) ' dPUFreq',dPUFreq
126     WRITE(20,*) ' a2MomXY',a2MomXY
127     WRITE(20,*) ' a2MomZ',a2MomZ
128     WRITE(20,*) ' a2TempXY',a2TempXY
129     WRITE(20,*) ' a2TempZ',a2TempZ
130     C********* GMGS *******************
131     C Ratio of time step
132     asyncfac = 1.0
133     C Scale dependent mixing scheme parameters
134     difref = a2tempXY
135     difmin = 0.1 * a2tempXY
136     difhmx = 2.5 * a2tempXY
137     coef = 1.8E9
138     dbotm = 1000.
139     C NOTE: MOD( INT(dslpcalcfreq), INT(dkcalcfreq) ) must = 0
140     dslpcalcfreq = 2. * delT
141     dkcalcfreq = 12. * delT
142     C Threashold for mix layer base
143     MixLayrDensGrad = 0.05
144     C Slope factors
145     slpmax1=5.e-3
146     slpmax2=1.e-2
147     slpmax3=1.e-3
148     small=1.e-30
149     dthin=50.
150     C Minimum density stratification
151     MinSigStrat= -1.e-8
152     c epsl = Av/Ah in Redis scheme
153     epsl = 0.
154     trDifScheme = gmconKh
155     verticalProfKh = .FALSE.
156    
157     divHmax = 1.D+300
158     divHMin = 0.
159     rDelt = 1./DELT
160     numberOfTimeSteps = NINT((endTime-startTime)/delT)
161     WRITE(20,*) ' numberOfTimeSteps',numberOfTimeSteps
162     currentTime = startTime
163     nIter = NINT(startTime/delt)
164     WRITE(20,*) ' nIter',nIter
165     WRITE(puSuffix,'(I10.10)') nIter
166     PickupRun = .FALSE.
167     IF ( startTime .NE. 0. ) pickupRun = .TRUE.
168     C State variables
169     gsNM1 = 0.
170     gtNM1 = 0.
171     guNM1 = 0.
172     gvNM1 = 0.
173     U = 0.
174     V = 0.
175     W = 0.
176     S = 0.
177     PH = 0.
178     PS = 0.
179     PSNM1 = 0.
180     PHNM1 = 0.
181     uAja = 0.
182     vAja = 0.
183     CENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNORE
184    
185    
186     THS(1) = 20.D0
187     THS(2) = 10.D0
188     C THS(3) = 10.D0
189     C THS(4) = 6.D0
190     SSPPT = 10.D0
191     T(_I3(1,:,:)) = 20.0D0
192     T(_I3(2,:,:)) = 10.0D0
193     C T(_I3(3,:,:)) = 10.0D0
194     C T(_I3(4,:,:)) = 6.0D0
195     C T(_I3(1,:,:)) = 1.D-12
196     C T(_I3(2,:,:)) = 0.
197     C T(_I3(3,:,:)) = 0.
198     C T(_I3(4,:,:)) = 0.
199     C T(_I3(1,:,:)) = 23.0D0
200     C T(_I3(2,:,:)) = 22.0D0
201     C T(_I3(3,:,:)) = 21.0D0
202     C T(_I3(4,:,:)) = 20.0D0
203     S = 10.0D0
204     DELTAZ(1) = 500.0D0
205     DELTAZ(2) = 500.0D0
206     C DELTAZ(3) = 1500.0D0
207     C DELTAZ(4) = 2500.0D0
208    
209     C Geometry
210     DM = 20. D3
211     WRITE(20,*) ' DM',DM
212    
213    
214     CBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNORE
215     rPvolHat = 0.
216     DO K = 1, nz
217     DELPS(K) = G*RONIL*DELTAZ(K)
218     XA(_I3(K,:,:)) = DELPS(K)*DM
219     YA(_I3(K,:,:)) = DELPS(K)*DM
220     rUVol(_I3(K,:,:)) = 1.D0/DM/DM/DELPS(K)
221     rVVol(_I3(K,:,:)) = 1.D0/DM/DM/DELPS(K)
222     rPVol(_I3(K,:,:)) = 1.D0/DM/DM/DELPS(K)
223     rPvolHat = rPvolHat+DM*DM*DELPS(K)
224     IF ( K .EQ. 1 ) THEN
225     uDdzOp(_I3(K,:,:)) = 1.D0/DELPS(K)
226     vDdzOp(_I3(K,:,:)) = 1.D0/DELPS(K)
227     pDdzOp(_I3(K,:,:)) = 1.D0/DELPS(K)
228     ELSE
229     uDdzOp(_I3(K,:,:)) = 2.D0/(DELPS(K)+DELPS(K-1))
230     vDdzOp(_I3(K,:,:)) = 2.D0/(DELPS(K)+DELPS(K-1))
231     pDdzOp(_I3(K,:,:)) = 2.D0/(DELPS(K)+DELPS(K-1))
232     ENDIF
233     ! begin *** Redi Operators
234     rDZatP(K)=1.D0/DELPS(K)
235     if (K.eq.1) then
236     rDZatW(K)=1.D0/(DELPS(K))
237     else
238     rDZatW(K)=2.D0/(DELPS(K)+DELPS(K-1))
239     endif
240     ! end ***** Redi Operators
241     ENDDO
242     ! begin *** Redi Operators
243     rDXatU=1.d0/DM
244     rDYatV=1.d0/DM
245     ! end ***** Redi Operators
246    
247     ZA = DM*DM
248     rPVolHat = 1.D0/rPvolHat
249     uDdxOp = 1.D0/DM
250     uDdyOp = 1.D0/DM
251     vDdxOp = uDdxOp
252     vDdyOp = uDdyOp
253     pDdxOp = uDdxOp
254     pDdyOp = uDdyOp
255     C Physical parameters
256     uTphiOp = 0.
257     vTphiOp = 0.
258     fCorW = 0.
259     pBeta = 1.D-11
260     CENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNOREENDIGNORE
261    
262    
263     C radial BETA code begins - to use comment out 'NOBETA CODE' following
264     XCENTER = (Nx + 1)/2
265     WRITE(0,*) ' XCENTER: ', XCENTER
266     YCENTER = (Ny + 1)/2
267     WRITE(0,*) ' YCENTER: ', YCENTER
268     RMAX = 54.0
269     WRITE(20,*) ' RMAX: ',RMAX
270     WRITE(20,*) ' actual radius (m): ',(RMAX*DM)
271     fRef = 4.0
272     WRITE(20,*) ' fRef: ',fRef
273     fWall = 4.0
274     WRITE(20,*) ' fWall: ',fWall
275     Beta = (fWall - fRef)/(RMAX)
276     WRITE(20,*) ' Beta(s-1/gridpoint): ',Beta
277     WRITE(20,*) ' Beta(s-1/m): ',Beta/DM
278    
279    
280    
281     C start NOBETACODE
282    
283     fCorV(_I3(:,:,1)) = 1.D-4
284     Cd2
285     C pBeta = 0.
286     Cd2
287     write(20,*) 'NO RADIAL BETA'
288    
289     fCorU(_I3(:,:,1)) = fCorV(_I3(:,:,1))+0.5*DM*pBeta
290     DO J = 2, ny
291     fCorV(_I3(:,:,J)) = fCorV(_I3(:,:,J-1))+DM*pBeta
292     fCorU(_I3(:,:,J)) = fCorU(_I3(:,:,J-1))+DM*pBeta
293     ENDDO
294    
295     C end NOBETACODE
296    
297    
298     CBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNOREBEGINIGNORE
299     sBeta = 0.
300     WRITE(20,*) ' sBeta',sBeta
301     tAlpha = 2.D-4
302     C tAlpha = 0.D-4
303     WRITE(20,*) ' tAlpha',tAlpha
304    
305    
306     C Numerics parameters
307     abEps = 0.1
308     rAja = 0.0D0
309     dUnit = 10
310     scrUnit1 = 11
311     scrUnit2 = 12
312     scrUnitError = 13
313     InitialisationError = .FALSE.
314     IntegrationError = .FALSE.
315    
316     C Topography
317     PMASK = WATER
318     C Enclose in a box
319     C Wall at X = Lx
320     PMASK(_I3(:,nx,:)) = LAND ! Land in east most cell
321     C Wall at Y = Ly
322     PMASK(_I3(:,:,ny)) = LAND ! Land in northern most cell
323     C bathySource = textBathymetry
324     C Island
325     C PMASK(_I3(:,1,24)) = LAND
326    
327     C
328    
329     MAX2DIT = 1000
330     TOLER2D = 1.D-13
331     freqCheckToler2d = 1
332    
333     C-- Momentum equation forcing terms
334     tauMax = 0.1D0
335     lY = 0.
336     DO j=1,nY-1
337     lY = lY + DM
338     ENDDO
339     distY = 0.
340     DO j=1,nY
341     distY = distY+DM*0.5D0
342     DO i=1,nX
343     fv(i,j) = 0.D0
344     fu(i,j) = -tauMax*cos(2.D0*PI*distY/lY)
345     C fu(i,j) = -tauMax
346     fu(i,j) = tauMax*sin(PI*distY/lY)
347     fu(i,j) = fu(i,j)/(deltaZ(1)*ronil)
348     ENDDO
349     distY = distY+DM*0.5D0
350     ENDDO
351     fu(4,4) = fu(4,4)*0.917D0
352    
353     C
354     CcnhDebugStarts
355     Cdbg WRITE(0,*) ' lY = ', lY
356     Cdbg WRITE(0,*) ' distY = ', distY
357     Cdbg WRITE(0,*) ' tauMax = ', tauMax
358     Cdbg WRITE(0,*) ' fu '
359     Cdbg CALL PLOT_FIELD(fu,Nx,Ny)
360     Cdbg STOP
361     CcnhDebugEnds
362     C
363     END

  ViewVC Help
Powered by ViewVC 1.1.22