/[MITgcm]/MITgcm/pkg/cheapaml/cheapaml_fields_load.F
ViewVC logotype

Annotation of /MITgcm/pkg/cheapaml/cheapaml_fields_load.F

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

Revision 1.11 - (hide annotations) (download)
Thu Apr 21 15:15:54 2011 UTC (13 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62w, checkpoint62y, checkpoint62x
Changes since 1.10: +69 -55 lines 
call S/R GET_PERIODIC_INTERVAL to get interp. weights and time record number
1 jmc 1.11 C $Header: /u/gcmpack/MITgcm/pkg/cheapaml/cheapaml_fields_load.F,v 1.10 2011/03/18 22:36:58 wienders Exp $
2 jmc 1.1 C $Name: $
3 jmc 1.7
4 jmc 1.1 #include "CHEAPAML_OPTIONS.h"
5 jmc 1.2
6 wienders 1.8
7 jmc 1.1 C !ROUTINE: CHEAPAML_FIELDS_LOAD
8     C !INTERFACE:
9     SUBROUTINE CHEAPAML_FIELDS_LOAD( myTime, myIter, myThid )
10     C *==========================================================*
11 jmc 1.2 C | SUBROUTINE CHEAPAML_FIELDS_LOAD
12     C | o Control reading of fields from external source.
13 jmc 1.1 C *==========================================================*
14    
15     C !USES:
16     IMPLICIT NONE
17     C === Global variables ===
18     #include "SIZE.h"
19     #include "EEPARAMS.h"
20     #include "PARAMS.h"
21     #include "FFIELDS.h"
22     c #include "GRID.h"
23     c #include "DYNVARS.h"
24     C #include "BULKF.h"
25 wienders 1.8 c#ifdef ALLOW_THSICE
26     c#include "THSICE_VARS.h"
27     c#endif
28 jmc 1.1 #include "CHEAPAML.h"
29 jmc 1.2
30 jmc 1.1 C !INPUT/OUTPUT PARAMETERS:
31     C === Routine arguments ===
32     C myThid - Thread no. that called this routine.
33     C myTime - Simulation time
34     C myIter - Simulation timestep number
35     C dsolms - Solar variation at Southern boundary
36     C dsolmn - Solar variation at Northern boundary
37     c xphaseinit - user input initial phase of year relative
38     c to mid winter. E.G. xphaseinit = pi implies time zero
39 jmc 1.2 c is mid summer.
40 jmc 1.1 INTEGER myThid
41     _RL myTime
42 jmc 1.7 _RL local
43 jmc 1.1 c _RL dsolms,dsolmn
44     c _RL xphaseinit
45     INTEGER myIter
46 jmc 1.2
47 jmc 1.1 C !LOCAL VARIABLES:
48     C === Local arrays ===
49     C trair[01] :: Relaxation temp. profile for air temperature
50     C qrair[01] :: Relaxation specific humidity profile for air
51     C solar[01] :: short wave flux
52     C uwind[01] :: zonal wind
53     C vwind[01] :: meridional wind
54 jmc 1.11 C CheaptracerR[01] :: Relaxation profile for passive tracer
55 jmc 1.7 C aWght, bWght :: Interpolation weights
56 jmc 1.1
57 jmc 1.11 COMMON /CHEAPAML_FIELDS_LOAD_I/ cheapaml_ldRec
58     COMMON /CHEAPAML_FIELDS_LOAD_RL/
59 jmc 1.7 & trair0, trair1,
60     & qrair0, qrair1,
61     & Solar0, Solar1,
62     & uwind0, uwind1,
63     & vwind0, vwind1,
64     & ustress0, ustress1,
65     & vstress0, vstress1,
66     & wavesh0, wavesh1,
67     & wavesp0, wavesp1,
68 wienders 1.9 & rair0, rair1,
69 jmc 1.11 & CheaptracerR0, CheaptracerR1
70 wienders 1.8
71 jmc 1.11 INTEGER cheapaml_ldRec(nSx,nSy)
72 jmc 1.7 _RL trair0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
73     _RL trair1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
74     _RL qrair0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
75     _RL qrair1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76     _RL Solar0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
77     _RL Solar1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
78     _RL uwind0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
79     _RL uwind1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
80     _RL vwind0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
81     _RL vwind1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
82     _RL ustress0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
83     _RL ustress1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
84     _RL vstress0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
85     _RL vstress1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
86     _RL wavesh0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
87     _RL wavesh1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
88     _RL wavesp0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
89     _RL wavesp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
90     _RL rair0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
91     _RL rair1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
92 wienders 1.9 _RL CheaptracerR0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
93     _RL CheaptracerR1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
94 jmc 1.1
95 jmc 1.11 INTEGER intimeP, intime0, intime1
96     INTEGER bi,bj,i,j
97     INTEGER iG,jG
98     _RL aWght,bWght,u
99 wienders 1.8 _RL ssq0,ssq1,ssq2,ssqa
100 jmc 1.2 c xsolph - phase of year, assuming time zero is mid winter
101 jmc 1.1 c xinxx - cos ( xsolph )
102 wienders 1.8 _RL xsolph,xinxx
103 jmc 1.1 c coefficients used to compute saturation specific humidity
104 wienders 1.8 DATA ssq0, ssq1, ssq2
105 jmc 1.1 & / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
106    
107 wienders 1.8 IF ( periodicExternalForcing_cheap ) THEN
108    
109 jmc 1.1
110     c the objective here is to give cheapaml a default periodic forcing
111     c consisting only of annually varying solar forcing, and thus Trelaxation
112     c variation. everything else, relative humidity, wind, are fixed. This
113 jmc 1.2 c keys off of solardata. if a solar data file exists, the model will
114 jmc 1.1 c assume there are files to be read and interpolated between, as is standard
115     c for the MITGCM.
116    
117 wienders 1.8 IF ( SolarFile .EQ. ' ' ) THEN
118     IF (useStressOption)then
119     write(*,*)' stress option is turned on. this is not
120     & consistent with the default time dependent forcing option'
121     stop
122 jmc 1.7 ENDIF
123 wienders 1.8 if ( myIter .EQ. nIter0 )then
124 jmc 1.1 WRITE(*,*)
125 wienders 1.8 & 'S/R Assuming Standard Annually Varying Solar Forcing'
126     endif
127     xsolph=myTime*2.d0*3.14159 _d 0/365. _d 0/86400. _d 0
128     xinxx=cos(xsolph+xphaseinit+3.14159 _d 0)
129     DO bj = myByLo(myThid), myByHi(myThid)
130     DO bi = myBxLo(myThid), myBxHi(myThid)
131     DO j=1,sNy
132     DO i=1,sNx
133     jG = myYGlobalLo-1+(bj-1)*sNy+j
134 jmc 1.1 local=225.d0+dsolms*xinxx-float((jg-1))/float((ny-1))*
135     & (37.5d0-dsolmn*xinxx)
136 wienders 1.8 Solar(i,j,bi,bj) = local
137     ENDDO
138 jmc 1.3 ENDDO
139     ENDDO
140 jmc 1.1 ENDDO
141 wienders 1.9 _EXCH_XY_RL(solar, mythid)
142 jmc 1.1 c relaxation temperature in radiative equilibrium
143 wienders 1.8 DO bj = myByLo(myThid), myByHi(myThid)
144     DO bi = myBxLo(myThid), myBxHi(myThid)
145     DO j=1,sNy
146     DO i=1,sNx
147     jG = myYGlobalLo-1+(bj-1)*sNy+j
148     local=solar(i,j,bi,bj)
149     local=(2.d0*local/stefan)**(0.25d0)-Celsius2K
150     TR(i,j,bi,bj) = local
151     ENDDO
152 jmc 1.3 ENDDO
153     ENDDO
154 jmc 1.1 ENDDO
155 wienders 1.9 _EXCH_XY_RL(TR, mythid)
156 jmc 1.1 c default specific humidity profile to 80% relative humidity
157 wienders 1.8 DO bj = myByLo(myThid), myByHi(myThid)
158     DO bi = myBxLo(myThid), myBxHi(myThid)
159     DO j=1,sNy
160     DO i=1,sNx
161     c jG = myYGlobalLo-1+(bj-1)*sNy+j
162     local = Tr(i,j,bi,bj)+Celsius2K
163     ssqa = ssq0*exp( lath*(ssq1-ssq2/local)) / p0
164     qr(i,j,bi,bj) = 0.8d0*ssqa
165     ENDDO
166 jmc 1.3 ENDDO
167     ENDDO
168 jmc 1.1 ENDDO
169 wienders 1.9 _EXCH_XY_RL(qr, mythid)
170 jmc 1.1 c u wind field
171 wienders 1.8 DO bj = myByLo(myThid), myByHi(myThid)
172     DO bi = myBxLo(myThid), myBxHi(myThid)
173     DO j=1,sNy
174     DO i=1,sNx
175     jG = myYGlobalLo-1+(bj-1)*sNy+j
176 jmc 1.11 local=-5.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))
177 wienders 1.8 uwind(i,j,bi,bj) = local
178     ENDDO
179 jmc 1.3 ENDDO
180     ENDDO
181 jmc 1.1 ENDDO
182 wienders 1.9 _EXCH_XY_RL(uwind, mythid)
183 jmc 1.1 c v wind field
184 wienders 1.8 DO bj = myByLo(myThid), myByHi(myThid)
185     DO bi = myBxLo(myThid), myBxHi(myThid)
186     DO j=1,sNy
187     DO i=1,sNx
188     jG = myYGlobalLo-1+(bj-1)*sNy+j
189     vwind(i,j,bi,bj) = 0.d0
190     ENDDO
191 jmc 1.7 ENDDO
192     ENDDO
193     ENDDO
194 wienders 1.9 _EXCH_XY_RL(vwind, mythid)
195     C Tracer field
196     IF (useCheapTracer) THEN
197     DO bj = myByLo(myThid), myByHi(myThid)
198     DO bi = myBxLo(myThid), myBxHi(myThid)
199     DO j=1,sNy
200     DO i=1,sNx
201     jG = myYGlobalLo-1+(bj-1)*sNy+j
202     CheaptracerR(i,j,bi,bj) = 290. _d 0
203     ENDDO
204     ENDDO
205     ENDDO
206     ENDDO
207     ENDIF
208 jmc 1.11 _EXCH_XY_RL(CheaptracerR, mythid)
209 wienders 1.9
210 wienders 1.8 ELSE
211 jmc 1.7
212 wienders 1.8 c here for usual interpolative forcings
213 jmc 1.7 C First call requires that we initialize everything to zero for safety
214 wienders 1.8 IF ( myIter .EQ. nIter0 ) THEN
215     DO bj = myByLo(myThid), myByHi(myThid)
216     DO bi = myBxLo(myThid), myBxHi(myThid)
217 jmc 1.11 cheapaml_ldRec(bi,bj) = 0
218 wienders 1.8 DO j=1-OLy,sNy+OLy
219     DO i=1-OLx,sNx+OLx
220     trair0 (i,j,bi,bj) = 0.
221     trair1 (i,j,bi,bj) = 0.
222     qrair0 (i,j,bi,bj) = 0.
223     qrair1 (i,j,bi,bj) = 0.
224     solar0 (i,j,bi,bj) = 0.
225     solar1 (i,j,bi,bj) = 0.
226     uwind0 (i,j,bi,bj) = 0.
227     uwind1 (i,j,bi,bj) = 0.
228     vwind0 (i,j,bi,bj) = 0.
229     vwind1 (i,j,bi,bj) = 0.
230     ustress0(i,j,bi,bj) = 0.
231     ustress1(i,j,bi,bj) = 0.
232     vstress0(i,j,bi,bj) = 0.
233     vstress1(i,j,bi,bj) = 0.
234     wavesh0 (i,j,bi,bj) = 0.
235     wavesh1 (i,j,bi,bj) = 0.
236     wavesp0 (i,j,bi,bj) = 0.
237     wavesp1 (i,j,bi,bj) = 0.
238     rair0 (i,j,bi,bj) = 0.
239     rair1 (i,j,bi,bj) = 0.
240 wienders 1.9 CheaptracerR0 (i,j,bi,bj) = 0.
241     CheaptracerR1 (i,j,bi,bj) = 0.
242 jmc 1.3 ENDDO
243     ENDDO
244 jmc 1.1 ENDDO
245 jmc 1.11 ENDDO
246 wienders 1.8 ENDIF
247 jmc 1.1
248 jmc 1.11 C-- Now calculate whether it is time to update the forcing arrays
249     CALL GET_PERIODIC_INTERVAL(
250     O intimeP, intime0, intime1, bWght, aWght,
251     I externForcingCycle_cheap, externForcingPeriod_cheap,
252     I deltaTclock, myTime, myThid )
253    
254     bi = myBxLo(myThid)
255     bj = myByLo(myThid)
256     #ifdef ALLOW_DEBUG
257     IF ( debugLevel.GE.debLevA ) THEN
258     _BEGIN_MASTER(myThid)
259     WRITE(standardMessageUnit,'(A,I10,A,4I5,A,2F14.10)')
260     & ' CHEAPAML_FIELDS_LOAD,', myIter,
261     & ' : iP,iLd,i0,i1=', intimeP,cheapaml_ldRec(bi,bj),
262     & intime0,intime1, ' ; Wght=', bWght, aWght
263     _END_MASTER(myThid)
264     ENDIF
265     #endif /* ALLOW_DEBUG */
266    
267     C- Make no assumption on sequence of calls to CHEAPAML_FIELDS_LOAD ;
268     C This is the correct formulation (works in Adjoint run).
269     IF ( intime1.NE.cheapaml_ldRec(bi,bj) ) THEN
270    
271     C-- If the above condition is met then we need to read in
272     C data for the period ahead and the period behind myTime.
273     _BEGIN_MASTER(myThid)
274     WRITE(standardMessageUnit,'(A,I10,A,2(2I5,A))')
275     & ' CHEAPAML_FIELDS_LOAD, it=', myIter,
276     & ' : Reading new data, i0,i1=', intime0, intime1,
277     & ' (prev=', intimeP, cheapaml_ldRec(bi,bj), ' )'
278     _END_MASTER(myThid)
279    
280 wienders 1.8 IF ( SolarFile .NE. ' ' ) THEN
281     CALL READ_REC_XY_RL( SolarFile,solar0,intime0,
282     & myIter,myThid )
283     CALL READ_REC_XY_RL( SolarFile,solar1,intime1,
284     & myIter,myThid )
285     ENDIF
286     IF ( TrFile .NE. ' ' ) THEN
287     CALL READ_REC_XY_RL( TRFile,trair0,intime0,
288     & myIter,myThid )
289     CALL READ_REC_XY_RL( TRFile,trair1,intime1,
290     & myIter,myThid )
291     ENDIF
292     IF ( QrFile .NE. ' ' ) THEN
293     CALL READ_REC_XY_RL( QrFile,qrair0,intime0,
294     & myIter,myThid )
295     CALL READ_REC_XY_RL( QrFile,qrair1,intime1,
296     & myIter,myThid )
297     ENDIF
298     IF ( UWindFile .NE. ' ' ) THEN
299     CALL READ_REC_XY_RL( UWindFile,uwind0,intime0,
300     & myIter,myThid )
301     CALL READ_REC_XY_RL( UWindFile,uwind1,intime1,
302     & myIter,myThid )
303     ENDIF
304     IF ( VWindFile .NE. ' ' ) THEN
305     CALL READ_REC_XY_RL( VWindFile,vwind0,intime0,
306     & myIter,myThid )
307     CALL READ_REC_XY_RL( VWindFile,vwind1,intime1,
308     & myIter,myThid )
309     ENDIF
310     IF(useStressOption)THEN
311     IF ( UStressFile .NE. ' ' ) THEN
312     CALL READ_REC_XY_RL( UStressFile,ustress0,intime0,
313     & myIter,myThid )
314     CALL READ_REC_XY_RL( UStressFile,ustress1,intime1,
315     & myIter,myThid )
316     ENDIF
317     IF ( VStressFile .NE. ' ' ) THEN
318     CALL READ_REC_XY_RL( VStressFile,vstress0,intime0,
319     & myIter,myThid )
320     CALL READ_REC_XY_RL( VStressFile,vstress1,intime1,
321     & myIter,myThid )
322     ENDIF
323     ENDIF
324     IF ( FluxFormula.eq.'COARE3') THEN
325     IF ( WaveHFile .NE. ' ' ) THEN
326     CALL READ_REC_XY_RL( WaveHFile,wavesh0,intime0,
327     & myIter,myThid )
328     CALL READ_REC_XY_RL( WaveHFile,wavesh1,intime1,
329     & myIter,myThid )
330     ENDIF
331     IF ( WavePFile .NE. ' ' ) THEN
332     CALL READ_REC_XY_RL( WavePFile,wavesp0,intime0,
333     & myIter,myThid )
334     CALL READ_REC_XY_RL( WavePFile,wavesp1,intime1,
335     & myIter,myThid )
336     ENDIF
337     ENDIF
338 wienders 1.9 IF (useCheapTracer) THEN
339     IF ( TracerRFile .NE. ' ' ) THEN
340     CALL READ_REC_XY_RL( TracerRFile,CheaptracerR0,intime0,
341     & myIter,myThid )
342     CALL READ_REC_XY_RL( TracerRFile,CheaptracerR1,intime1,
343     & myIter,myThid )
344 jmc 1.11 ENDIF
345 wienders 1.9 ENDIF
346 wienders 1.8 _EXCH_XY_RL(trair0 , myThid )
347     _EXCH_XY_RL(qrair0 , myThid )
348     _EXCH_XY_RL(solar0 , myThid )
349     _EXCH_XY_RL(uwind0 , myThid )
350     _EXCH_XY_RL(vwind0 , myThid )
351     _EXCH_XY_RL(trair1 , myThid )
352     _EXCH_XY_RL(qrair1 , myThid )
353     _EXCH_XY_RL(solar1 , myThid )
354     _EXCH_XY_RL(uwind1 , myThid )
355     _EXCH_XY_RL(vwind1 , myThid )
356     IF(useStressOption)THEN
357     _EXCH_XY_RL(uwind0 , myThid )
358     _EXCH_XY_RL(vwind0 , myThid )
359     _EXCH_XY_RL(uwind1 , myThid )
360     _EXCH_XY_RL(vwind1 , myThid )
361     ENDIF
362     IF(FluxFormula.EQ.'COARE3') THEN
363 wienders 1.9 _EXCH_XY_RL(wavesp0 , myThid )
364     _EXCH_XY_RL(wavesp1 , myThid )
365     _EXCH_XY_RL(wavesh0 , myThid )
366     _EXCH_XY_RL(wavesh1 , myThid )
367     ENDIF
368     IF(useCheapTracer) THEN
369     _EXCH_XY_RL(CheaptracerR0 , myThid )
370     _EXCH_XY_RL(CheaptracerR1, myThid )
371 jmc 1.11 ENDIF
372    
373     C- save newly loaded time-record
374     DO bj = myByLo(myThid), myByHi(myThid)
375     DO bi = myBxLo(myThid), myBxHi(myThid)
376     cheapaml_ldRec(bi,bj) = intime1
377     ENDDO
378     ENDDO
379    
380 wienders 1.8 C end of loading new fields block
381     ENDIF
382 jmc 1.1
383 wienders 1.8 C-- Interpolate TR, QR, SOLAR
384     DO bj = myByLo(myThid), myByHi(myThid)
385     DO bi = myBxLo(myThid), myBxHi(myThid)
386     DO j=1-OLy,sNy+OLy
387     DO i=1-OLx,sNx+OLx
388     TR(i,j,bi,bj) = bWght*trair0(i,j,bi,bj)
389     & +aWght*trair1(i,j,bi,bj) !+273.15
390     qr(i,j,bi,bj) = bWght*qrair0(i,j,bi,bj)
391     & +aWght*qrair1(i,j,bi,bj)
392     uwind(i,j,bi,bj) = bWght*uwind0(i,j,bi,bj)
393     & +aWght*uwind1(i,j,bi,bj)
394     vwind(i,j,bi,bj) = bWght*vwind0(i,j,bi,bj)
395     & +aWght*vwind1(i,j,bi,bj)
396     solar(i,j,bi,bj) = bWght*solar0(i,j,bi,bj)
397     & +aWght*solar1(i,j,bi,bj)
398     IF(useStressOption)THEN
399     ustress(i,j,bi,bj) = bWght*ustress0(i,j,bi,bj)
400 jmc 1.7 & +aWght*ustress1(i,j,bi,bj)
401 wienders 1.8 vstress(i,j,bi,bj) = bWght*vstress0(i,j,bi,bj)
402 jmc 1.7 & +aWght*vstress1(i,j,bi,bj)
403 wienders 1.8 ENDIF
404 wienders 1.9 IF(useCheapTracer)THEN
405     CheaptracerR(i,j,bi,bj) = bWght*CheaptracerR0(i,j,bi,bj)
406     & +aWght*CheaptracerR1(i,j,bi,bj)
407 jmc 1.11 ENDIF
408 wienders 1.8 IF(FluxFormula.eq.'COARE3')THEN
409     IF(WaveHFile.ne.' ')THEN
410     wavesh(i,j,bi,bj) = bWght*wavesh0(i,j,bi,bj)
411     & +aWght*wavesh1(i,j,bi,bj)
412     ENDIF
413     IF(WavePFile.ne.' ')THEN
414     wavesp(i,j,bi,bj) = bWght*wavesp0(i,j,bi,bj)
415     & +aWght*wavesp1(i,j,bi,bj)
416     ENDIF
417     ELSE
418     u=uwind(i,j,bi,bj)**2+vwind(i,j,bi,bj)**2
419     u=dsqrt(u)
420 jmc 1.11 wavesp(i,j,bi,bj)=0.729 _d 0 * u
421 wienders 1.8 wavesh(i,j,bi,bj)=0.018 _d 0 * u*u*(1. _d 0 + .015 _d 0 *u)
422     ENDIF
423 jmc 1.1 ENDDO
424     ENDDO
425 wienders 1.8 ENDDO
426     ENDDO
427    
428 jmc 1.11 c end if solarFile is empty
429 wienders 1.8 ENDIF
430 jmc 1.1
431 wienders 1.8 c end of periodic forcing options, on to steady option
432 jmc 1.7
433 jmc 1.1 ELSE
434 jmc 1.7
435 jmc 1.1 IF ( myIter .EQ. nIter0 ) THEN
436     IF ( SolarFile .NE. ' ' ) THEN
437 wienders 1.8 CALL READ_FLD_XY_RS( SolarFile,' ',solar,0,myThid )
438 jmc 1.3 ELSE
439     DO bj = myByLo(myThid), myByHi(myThid)
440     DO bi = myBxLo(myThid), myBxHi(myThid)
441     DO j=1,sNy
442     DO i=1,sNx
443 jmc 1.1 jG = myYGlobalLo-1+(bj-1)*sNy+j
444     local=225.d0-float((jg-1))/float((ny-1))*37.5d0
445     Solar(i,j,bi,bj) = local
446     ENDDO
447 jmc 1.3 ENDDO
448     ENDDO
449 jmc 1.1 ENDDO
450 jmc 1.3 ENDIF
451 wienders 1.8 _EXCH_XY_RS(solar, mythid)
452 jmc 1.1 IF ( TrFile .NE. ' ' ) THEN
453 jmc 1.7 CALL READ_FLD_XY_RL( TrFile,' ',tr,0,myThid )
454 jmc 1.3 ELSE
455     DO bj = myByLo(myThid), myByHi(myThid)
456     DO bi = myBxLo(myThid), myBxHi(myThid)
457     DO j=1,sNy
458     DO i=1,sNx
459 jmc 1.7 jG = myYGlobalLo-1+(bj-1)*sNy+j
460     local=solar(i,j,bi,bj)
461     local=(2.d0*local/stefan)**(0.25d0)-273.16
462 wienders 1.8 TR(i,j,bi,bj) = local
463 jmc 1.1 ENDDO
464 jmc 1.3 ENDDO
465     ENDDO
466 jmc 1.1 ENDDO
467     ENDIF
468 jmc 1.11
469 wienders 1.8 _EXCH_XY_RL(TR, mythid)
470     c do specific humidity
471     IF ( QrFile .NE. ' ') THEN
472 jmc 1.11 CALL READ_FLD_XY_RL( QrFile,' ',qr,0,myThid )
473 jmc 1.3 ELSE
474 wienders 1.8 c default specific humidity profile to 80% relative humidity
475 jmc 1.3 DO bj = myByLo(myThid), myByHi(myThid)
476     DO bi = myBxLo(myThid), myBxHi(myThid)
477     DO j=1,sNy
478     DO i=1,sNx
479 jmc 1.1 c jG = myYGlobalLo-1+(bj-1)*sNy+j
480     local = Tr(i,j,bi,bj)+273.16d0
481     ssqa = ssq0*exp( lath*(ssq1-ssq2/local)) / p0
482     qr(i,j,bi,bj) = 0.8d0*ssqa
483     ENDDO
484 jmc 1.3 ENDDO
485     ENDDO
486 jmc 1.1 ENDDO
487     ENDIF
488 wienders 1.8 _EXCH_XY_RL(qr, mythid)
489 jmc 1.1 IF ( UWindFile .NE. ' ' ) THEN
490 jmc 1.7 CALL READ_FLD_XY_RL( UWindFile,' ',uwind,0,myThid )
491 jmc 1.3 ELSE
492     DO bj = myByLo(myThid), myByHi(myThid)
493     DO bi = myBxLo(myThid), myBxHi(myThid)
494     DO j=1,sNy
495     DO i=1,sNx
496 wienders 1.8 jG = myYGlobalLo-1+(bj-1)*sNy+j
497 jmc 1.1 c mod for debug
498     c to return to original code, uncomment following line
499     c comment out 2nd line
500 wienders 1.8 local=-5.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))
501     c local=0.d0*cos(2.d0*pi*float(jg-1)/(float(ny-1)))
502     uwind(i,j,bi,bj) = local
503 jmc 1.1 ENDDO
504 jmc 1.3 ENDDO
505     ENDDO
506 jmc 1.1 ENDDO
507     ENDIF
508 wienders 1.8 _EXCH_XY_RL(uwind, mythid)
509 jmc 1.1 IF ( VWindFile .NE. ' ' ) THEN
510 jmc 1.7 CALL READ_FLD_XY_RL( VWindFile,' ',vwind,0,myThid )
511 jmc 1.3 ELSE
512     DO bj = myByLo(myThid), myByHi(myThid)
513     DO bi = myBxLo(myThid), myBxHi(myThid)
514     DO j=1,sNy
515     DO i=1,sNx
516 wienders 1.8 jG = myYGlobalLo-1+(bj-1)*sNy+j
517     vwind(i,j,bi,bj) = 0.d0
518 jmc 1.1 ENDDO
519 jmc 1.3 ENDDO
520     ENDDO
521 jmc 1.1 ENDDO
522     ENDIF
523 wienders 1.8 _EXCH_XY_RL(vwind, mythid)
524     IF(useStressOption)THEN
525     IF ( UStressFile .NE. ' ' ) THEN
526     CALL READ_FLD_XY_RL( UStressFile,' ',ustress,0,myThid )
527     ELSE
528     write(*,*)' U Stress File absent with stress option'
529     stop
530     ENDIF
531     IF ( VStressFile .NE. ' ' ) THEN
532 wienders 1.9 CALL READ_FLD_XY_RL( VStressFile,' ',vstress,0,myThid )
533 wienders 1.8 ELSE
534     write(*,*)' V Stress File absent with stress option'
535     stop
536     ENDIF
537     _EXCH_XY_RL(ustress, mythid)
538     _EXCH_XY_RL(vstress, mythid)
539     ENDIF
540 wienders 1.9 IF(useCheapTracer)THEN
541     IF ( TracerRFile .NE. ' ' ) THEN
542     CALL READ_FLD_XY_RL( TracerRFile,' ',CheaptracerR,0,myThid )
543     ELSE
544     DO bj = myByLo(myThid), myByHi(myThid)
545     DO bi = myBxLo(myThid), myBxHi(myThid)
546     DO j=1,sNy
547     DO i=1,sNx
548     CheaptracerR(i,j,bi,bj)=290. _d 0
549     ENDDO
550     ENDDO
551     ENDDO
552     ENDDO
553     ENDIF
554     _EXCH_XY_RL(CheaptracerR, mythid)
555 jmc 1.11 ENDIF
556 wienders 1.8 IF (FluxFormula.eq.'COARE3')THEN
557     IF (WaveHFile.NE.' ')THEN
558     CALL READ_FLD_XY_RL( WaveHFile,' ',wavesh,0,myThid )
559     ENDIF
560     IF (WavePFile.NE.' ')THEN
561     CALL READ_FLD_XY_RL( WavePFile,' ',wavesp,0,myThid )
562     ELSE
563     DO bj = myByLo(myThid), myByHi(myThid)
564     DO bi = myBxLo(myThid), myBxHi(myThid)
565     DO j=1,sNy
566     DO i=1,sNx
567     u=uwind(i,j,bi,bj)**2+vwind(i,j,bi,bj)**2
568     u=dsqrt(u)
569 jmc 1.11 wavesp(i,j,bi,bj)=0.729 _d 0 * u
570 wienders 1.8 wavesh(i,j,bi,bj)=0.018 _d 0 * u*u*(1. _d 0 + .015 _d 0 *u)
571 jmc 1.7 ENDDO
572     ENDDO
573     ENDDO
574 wienders 1.8 ENDDO
575     ENDIF
576     _EXCH_XY_RL(wavesp, mythid)
577     _EXCH_XY_RL(wavesh, mythid)
578 jmc 1.11 ENDIF
579 wienders 1.8 ENDIF
580 jmc 1.7
581 jmc 1.2
582 jmc 1.7 C endif for Steady Option
583 jmc 1.1 ENDIF
584    
585 jmc 1.7 C fill in outer edges
586 wienders 1.8
587     DO bj = myByLo(myThid), myByHi(myThid)
588     DO bi = myBxLo(myThid), myBxHi(myThid)
589     do j=1-oly,sny+oly
590     jG = myYGlobalLo-1+(bj-1)*sNy+j
591     do i=1-olx,snx+olx
592     iG=myXGlobalLo-1+(bi-1)*sNx+i
593     if(iG.lt.1)then
594     Tr(i,j,bi,bj)=Tr(1,j,bi,bj)
595     qr(i,j,bi,bj)=qr(1,j,bi,bj)
596     uwind(i,j,bi,bj)=uwind(1,j,bi,bj)
597     vwind(i,j,bi,bj)=vwind(1,j,bi,bj)
598     Solar(i,j,bi,bj)=Solar(1,j,bi,bj)
599 jmc 1.11
600 wienders 1.9 IF(useStressOption)THEN
601 wienders 1.8 ustress(i,j,bi,bj)=ustress(1,j,bi,bj)
602     vstress(i,j,bi,bj)=vstress(1,j,bi,bj)
603 wienders 1.9 ENDIF
604    
605     IF(useCheapTracer)THEN
606     CheaptracerR(i,j,bi,bj)=CheaptracerR(1,j,bi,bj)
607 jmc 1.11 ENDIF
608 wienders 1.8
609     if(FluxFormula.eq.'COARE3')then
610     wavesp(i,j,bi,bj)=wavesp(1,j,bi,bj)
611     wavesh(i,j,bi,bj)=wavesh(1,j,bi,bj)
612 jmc 1.7 endif
613 wienders 1.8
614     elseif(iG.gt.Nx)then
615     Tr(i,j,bi,bj)=Tr(sNx,j,bi,bj)
616     qr(i,j,bi,bj)=qr(sNx,j,bi,bj)
617     uwind(i,j,bi,bj)=uwind(sNx,j,bi,bj)
618     vwind(i,j,bi,bj)=vwind(sNx,j,bi,bj)
619     Solar(i,j,bi,bj)=Solar(sNx,j,bi,bj)
620    
621 jmc 1.7 if(UseStressOption)then
622 wienders 1.8 ustress(i,j,bi,bj)=ustress(sNx,j,bi,bj)
623     vstress(i,j,bi,bj)=vstress(sNx,j,bi,bj)
624 jmc 1.7 endif
625 wienders 1.8
626 wienders 1.9 IF(useCheapTracer)THEN
627     CheaptracerR(i,j,bi,bj)=CheaptracerR(sNx,j,bi,bj)
628 jmc 1.11 ENDIF
629 wienders 1.9
630 wienders 1.8 if(FluxFormula.eq.'COARE3')then
631     wavesp(i,j,bi,bj)=wavesp(sNx,j,bi,bj)
632     wavesh(i,j,bi,bj)=wavesh(sNx,j,bi,bj)
633 jmc 1.7 endif
634 wienders 1.8
635     elseif(jG.lt.1)then
636     Tr(i,j,bi,bj)=Tr(i,1,bi,bj)
637     qr(i,j,bi,bj)=qr(i,1,bi,bj)
638     uwind(i,j,bi,bj)=uwind(i,1,bi,bj)
639     vwind(i,j,bi,bj)=vwind(i,1,bi,bj)
640     Solar(i,j,bi,bj)=Solar(i,1,bi,bj)
641    
642 jmc 1.7 if(UseStressOption)then
643 wienders 1.8 ustress(i,j,bi,bj)=ustress(i,1,bi,bj)
644     vstress(i,j,bi,bj)=vstress(i,1,bi,bj)
645 jmc 1.7 endif
646 wienders 1.8
647 wienders 1.9 IF(useCheapTracer)THEN
648     CheaptracerR(i,j,bi,bj)=CheaptracerR(i,1,bi,bj)
649 jmc 1.11 ENDIF
650 wienders 1.9
651 wienders 1.8 if(FluxFormula.eq.'COARE3')then
652     wavesp(i,j,bi,bj)=wavesp(i,1,bi,bj)
653     wavesh(i,j,bi,bj)=wavesh(i,1,bi,bj)
654 jmc 1.7 endif
655 wienders 1.8
656     elseif(jG.gt.Ny)then
657     Tr(i,j,bi,bj)=Tr(i,sNy,bi,bj)
658     qr(i,j,bi,bj)=qr(i,sNy,bi,bj)
659     uwind(i,j,bi,bj)=uwind(i,sNy,bi,bj)
660     vwind(i,j,bi,bj)=vwind(i,sNy,bi,bj)
661     Solar(i,j,bi,bj)=Solar(i,sNy,bi,bj)
662    
663 jmc 1.7 if(UseStressOption)then
664 wienders 1.8 ustress(i,j,bi,bj)=ustress(i,sNy,bi,bj)
665     vstress(i,j,bi,bj)=vstress(i,sNy,bi,bj)
666 jmc 1.7 endif
667 wienders 1.8
668 wienders 1.9 IF(useCheapTracer)THEN
669     CheaptracerR(i,j,bi,bj)=CheaptracerR(i,sNy,bi,bj)
670 jmc 1.11 ENDIF
671 wienders 1.9
672 wienders 1.8 if(FluxFormula.eq.'COARE3')then
673     wavesp(i,j,bi,bj)=wavesp(i,sNy,bi,bj)
674     wavesh(i,j,bi,bj)=wavesh(i,sNy,bi,bj)
675 jmc 1.7 endif
676 wienders 1.8
677     endif
678 jmc 1.7 ENDDO
679 wienders 1.8 ENDDO
680 jmc 1.7 ENDDO
681 wienders 1.8 ENDDO
682     RETURN
683 jmc 1.1 END
  ViewVC Help
Powered by ViewVC 1.1.22