/[MITgcm]/MITgcm_contrib/ocean_inversion_project/code/ptracers_init.F
ViewVC logotype

Annotation of /MITgcm_contrib/ocean_inversion_project/code/ptracers_init.F

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


Revision 1.7 - (hide annotations) (download)
Wed Oct 22 16:02:04 2003 UTC (21 years, 9 months ago) by dimitri
Branch: MAIN
Changes since 1.6: +217 -30 lines
added netcdf support to ptracers_forcing_surf.F

1 dimitri 1.7 C $Header: /usr/local/gcmpack/MITgcm_contrib/ocean_inversion_project/code/ptracers_init.F,v 1.6 2003/10/21 07:31:11 dimitri Exp $
2 dimitri 1.1 C $Name: $
3    
4     #include "PTRACERS_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: PTRACERS_INIT
8    
9     C !INTERFACE: ==========================================================
10     SUBROUTINE PTRACERS_INIT( myThid )
11    
12     C !DESCRIPTION:
13     C Initialize PTRACERS data structures
14 dimitri 1.2 C This file is customized to compute CO2 perturbations from
15     C 30 ocean regions for Gruber's ocean inversion project.
16 dimitri 1.1
17     C !USES: ===============================================================
18     IMPLICIT NONE
19     #include "SIZE.h"
20     #include "EEPARAMS.h"
21     #include "PARAMS.h"
22     #include "GRID.h"
23     #include "PTRACERS.h"
24    
25    
26     C !INPUT PARAMETERS: ===================================================
27     C myThid :: thread number
28     INTEGER myThid
29    
30     C !OUTPUT PARAMETERS: ==================================================
31     C none
32    
33     #ifdef ALLOW_PTRACERS
34    
35     C !LOCAL VARIABLES: ====================================================
36     C i,j,k,bi,bj,iTracer :: loop indices
37 dimitri 1.4 INTEGER i,j,k,bi,bj,iTracer,iMonth,iUnit,iRec
38 dimitri 1.1 CHARACTER*(10) suff
39 dimitri 1.2 _RL SumPtracer
40 dimitri 1.5
41     #ifdef OCEAN_INVERSION_NETCDF
42 dimitri 1.7 Real*8 global(Nx,Ny)
43     _RL local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44 dimitri 1.5 Real lon(Nx,Ny),bounds_lon(Nx,Ny,2,2)
45     Real lat(Nx,Ny),bounds_lat(Nx,Ny,2,2)
46     Real depth(Nr),bounds_depth(Nr,2)
47     Real MASK(Nx,Ny,Nr)
48     Real AREA(Nx,Ny)
49     Real BATHY(Nx,Ny)
50     Real REGION_MASK(Nx,Ny,PTRACERS_num)
51     Real REGION_AREA(PTRACERS_num)
52     #endif /* OCEAN_INVERSION_NETCDF */
53    
54 dimitri 1.1 CEOP
55    
56     C Loop over tracers
57     DO iTracer = 1, PTRACERS_num
58    
59     C Loop over tiles
60 dimitri 1.2 DO bj = myByLo(myThid), myByHi(myThid)
61     DO bi = myBxLo(myThid), myBxHi(myThid)
62 dimitri 1.1
63     C Initialize arrays in common blocks :
64 dimitri 1.2 DO k=1,Nr
65     DO j=1-Oly,sNy+OLy
66     DO i=1-Olx,sNx+Olx
67     pTracer(i,j,k,bi,bj,iTracer) = 0. _d 0
68     gPtr(i,j,k,bi,bj,iTracer) = 0. _d 0
69     gPtrNM1(i,j,k,bi,bj,iTracer) = 0. _d 0
70     ENDDO
71 dimitri 1.3 ENDDO
72     ENDDO
73     DO j=1-Oly,sNy+OLy
74     DO i=1-Olx,sNx+Olx
75     surfaceTendencyPtr(i,j,bi,bj,iTracer) = 0. _d 0
76 dimitri 1.1 ENDDO
77     ENDDO
78    
79     C end bi,bj loops
80 dimitri 1.2 ENDDO
81 dimitri 1.1 ENDDO
82    
83     C end of Tracer loop
84     ENDDO
85    
86 dimitri 1.2 C Read from a pickup file if nIter0 is not zero
87     IF (nIter0.NE.0) THEN
88 dimitri 1.1 C-- Suffix for pickup files
89     IF (pickupSuff.EQ.' ') THEN
90     WRITE(suff,'(I10.10)') nIter0
91     ELSE
92     WRITE(suff,'(A10)') pickupSuff
93     ENDIF
94     CALL PTRACERS_READ_CHECKPOINT( nIter0,myThid )
95     ENDIF
96    
97 dimitri 1.2 C Initialize pTracerMasks
98 dimitri 1.1 CALL PTRACERS_READ_MASK ( mythid )
99 dimitri 1.2
100     C Initialize pTracerTakahashi
101     CALL PTRACERS_READ_Takahashi ( mythid )
102    
103     C Normalize pTracerTakahashi so that 1e18 mol/yr is released
104     C from each model region defined in pTracerMasks.
105     C It is assumed that each year is 365.25 days (31557600 s)
106     C long and that each month is 2629800 s.
107    
108     DO iTracer = 1, PTRACERS_num
109     SumPtracer = 0.
110     DO bj = myByLo(myThid), myByHi(myThid)
111     DO bi = myBxLo(myThid), myBxHi(myThid)
112     DO j=1,sNy
113     DO i=1,sNx
114     DO iMonth=1,12
115     SumPtracer = SumPtracer + 2629800. * rA(I,J,bi,bj) *
116     & pTracerMasks (I,J,iTracer,bi,bj) *
117     & pTracerTakahashi(I,J,iMonth ,bi,bj)
118     ENDDO
119     ENDDO
120     ENDDO
121     ENDDO
122     ENDDO
123     _GLOBAL_SUM_R8( SumPtracer, myThid )
124     DO bj = myByLo(myThid), myByHi(myThid)
125     DO bi = myBxLo(myThid), myBxHi(myThid)
126     DO j=1,sNy
127     DO i=1,sNx
128     DO iMonth=1,12
129     IF ( pTracerMasks(I,J,iTracer,bi,bj) .eq. 1. )
130     & pTracerTakahashi(I,J,iMonth ,bi,bj) = 1.e18 *
131     & pTracerTakahashi(I,J,iMonth ,bi,bj) / SumPtracer
132     ENDDO
133     ENDDO
134     ENDDO
135     ENDDO
136     ENDDO
137     ENDDO
138 dimitri 1.4
139     #ifdef OCEAN_INVERSION_PROJECT_TIME_DEPENDENT
140     C Initialize atmospheric CO2 array
141     call mdsfindunit( iUnit, mythid)
142     open(iUnit,file='splco2_cis92a.dat',status='old',form="FORMATTED")
143     do iRec=1,pTracerAtm_Nrec
144     read(iUnit,*) pTracerAtmYear(iRec), pTracerAtmCO2(iRec)
145     cdb print*, '###', iRec, pTracerAtmYear(iRec), pTracerAtmCO2(iRec)
146     enddo
147     close(iUnit)
148     #endif /* OCEAN_INVERSION_PROJECT_TIME_DEPENDENT */
149 dimitri 1.1
150 dimitri 1.5 #ifdef OCEAN_INVERSION_NETCDF
151 dimitri 1.7
152     C-- lon
153     DO bj = myByLo(myThid), myByHi(myThid)
154     DO bi = myBxLo(myThid), myBxHi(myThid)
155     DO j=1,sNy
156     DO i=1,sNx
157     local (i,j,bi,bj) = xC (i,j,bi,bj)
158     ENDDO
159     ENDDO
160     ENDDO
161     ENDDO
162     call gather_2d ( global, local, myThid )
163     DO j=1,Ny
164     DO i=1,Nx
165     if ( global(i,j) .gt. 180. ) then
166     lon (i,j) = global(i,j) - 360.
167     else
168     lon (i,j) = global(i,j)
169     endif
170     ENDDO
171     ENDDO
172    
173     C-- bounds_lon
174     DO bj = myByLo(myThid), myByHi(myThid)
175     DO bi = myBxLo(myThid), myBxHi(myThid)
176     DO j=1,sNy
177     DO i=1,sNx
178     local (i,j,bi,bj) = xG (i,j,bi,bj)
179     ENDDO
180     ENDDO
181     ENDDO
182     ENDDO
183     call gather_2d ( global, local, myThid )
184     DO j=1,Ny
185     DO i=1,Nx
186     if ( global(i,j) .gt. 180. ) then
187     bounds_lon (i,j,1,1) = global(i,j) - 360.
188     bounds_lon (i,j,1,2) = global(i,j) - 360.
189     else
190     bounds_lon (i,j,1,1) = global(i,j)
191     bounds_lon (i,j,1,2) = global(i,j)
192     endif
193     ENDDO
194     ENDDO
195     DO bj = myByLo(myThid), myByHi(myThid)
196     DO bi = myBxLo(myThid), myBxHi(myThid)
197     DO j=1,sNy
198     DO i=1,sNx
199     local (i,j,bi,bj) = xG (i+1,j,bi,bj)
200     ENDDO
201     ENDDO
202     ENDDO
203     ENDDO
204     call gather_2d ( global, local, myThid )
205     DO j=1,Ny
206     DO i=1,Nx
207     if ( global(i,j) .gt. 180. ) then
208     bounds_lon (i,j,2,1) = global(i,j) - 360.
209     bounds_lon (i,j,2,2) = global(i,j) - 360.
210     else
211     bounds_lon (i,j,2,1) = global(i,j)
212     bounds_lon (i,j,2,2) = global(i,j)
213     endif
214     ENDDO
215     ENDDO
216    
217     C-- lat
218     DO bj = myByLo(myThid), myByHi(myThid)
219     DO bi = myBxLo(myThid), myBxHi(myThid)
220     DO j=1,sNy
221     DO i=1,sNx
222     local (i,j,bi,bj) = yC (i,j,bi,bj)
223     ENDDO
224     ENDDO
225     ENDDO
226     ENDDO
227     call gather_2d ( global, local, myThid )
228     DO j=1,Ny
229     DO i=1,Nx
230     lat (i,j) = global(i,j)
231     ENDDO
232     ENDDO
233    
234     C-- bounds_lat
235     DO bj = myByLo(myThid), myByHi(myThid)
236     DO bi = myBxLo(myThid), myBxHi(myThid)
237     DO j=1,sNy
238     DO i=1,sNx
239     local (i,j,bi,bj) = yG (i,j,bi,bj)
240     ENDDO
241     ENDDO
242     ENDDO
243     ENDDO
244     call gather_2d ( global, local, myThid )
245     DO j=1,Ny
246     DO i=1,Nx
247     bounds_lat (i,j,1,1) = global(i,j)
248     bounds_lat (i,j,2,1) = global(i,j)
249     ENDDO
250     ENDDO
251     DO bj = myByLo(myThid), myByHi(myThid)
252     DO bi = myBxLo(myThid), myBxHi(myThid)
253     DO j=1,sNy
254     DO i=1,sNx
255     local (i,j,bi,bj) = yG (i,j+1,bi,bj)
256     ENDDO
257     ENDDO
258     ENDDO
259     ENDDO
260     call gather_2d ( global, local, myThid )
261     DO j=1,Ny
262     DO i=1,Nx
263     bounds_lat (i,j,1,2) = global(i,j)
264     bounds_lat (i,j,2,2) = global(i,j)
265     ENDDO
266     ENDDO
267    
268     C-- AREA
269     DO bj = myByLo(myThid), myByHi(myThid)
270     DO bi = myBxLo(myThid), myBxHi(myThid)
271     DO j=1,sNy
272     DO i=1,sNx
273     local (i,j,bi,bj) = rA (i,j,bi,bj)
274     ENDDO
275     ENDDO
276     ENDDO
277     ENDDO
278     call gather_2d ( global, local, myThid )
279 dimitri 1.5 DO j=1,Ny
280 dimitri 1.7 DO i=1,Nx
281     AREA (i,j) = global(i,j)
282     ENDDO
283     ENDDO
284    
285     C-- BATHY
286     DO bj = myByLo(myThid), myByHi(myThid)
287     DO bi = myBxLo(myThid), myBxHi(myThid)
288     DO j=1,sNy
289     DO i=1,sNx
290     local (i,j,bi,bj) = R_low (i,j,bi,bj)
291     ENDDO
292     ENDDO
293     ENDDO
294     ENDDO
295     call gather_2d ( global, local, myThid )
296     DO j=1,Ny
297     DO i=1,Nx
298     BATHY (i,j) = abs ( global(i,j) )
299     ENDDO
300     ENDDO
301    
302     C-- MASK
303     DO k=1,Nr
304     DO bj = myByLo(myThid), myByHi(myThid)
305     DO bi = myBxLo(myThid), myBxHi(myThid)
306     DO j=1,sNy
307     DO i=1,sNx
308     local (i,j,bi,bj) = hFacC (i,j,k,bi,bj)
309     ENDDO
310     ENDDO
311     ENDDO
312     ENDDO
313     call gather_2d ( global, local, myThid )
314     DO j=1,Ny
315     DO i=1,Nx
316     MASK (i,j,k) = global(i,j)
317     ENDDO
318     ENDDO
319     ENDDO
320    
321     C-- REGION_MASK
322     DO iTracer = 1, PTRACERS_num
323     DO bj = myByLo(myThid), myByHi(myThid)
324     DO bi = myBxLo(myThid), myBxHi(myThid)
325     DO j=1,sNy
326     DO i=1,sNx
327     local (i,j,bi,bj) = pTracerMasks(i,j,iTracer,bi,bj)
328     ENDDO
329 dimitri 1.5 ENDDO
330 dimitri 1.7 ENDDO
331     ENDDO
332     call gather_2d ( global, local, myThid )
333     DO j=1,Ny
334     DO i=1,Nx
335     REGION_MASK (i,j,iTracer) = global(i,j)
336     ENDDO
337     ENDDO
338 dimitri 1.5 ENDDO
339 dimitri 1.7
340     C-- depth, bounds_depth
341 dimitri 1.5 DO k=1,Nr
342 dimitri 1.7 depth (k ) = abs ( rC (k) )
343     bounds_depth (k,1) = abs ( rF (k+1) )
344     bounds_depth (k,2) = abs ( rF (k) )
345 dimitri 1.5 ENDDO
346 dimitri 1.7
347 dimitri 1.5 DO iTracer = 1, PTRACERS_num
348 dimitri 1.7 REGION_AREA (iTracer) = 0.
349     DO j=1,Ny
350     DO i=1,Nx
351     REGION_AREA (iTracer) = REGION_AREA (iTracer) +
352     & REGION_MASK (i,j,iTracer) * AREA (i,j)
353     ENDDO
354     ENDDO
355 dimitri 1.5 ENDDO
356 dimitri 1.7
357     C-- Master thread of process 0 writes netcdf output file
358     _BEGIN_MASTER( myThid )
359     #ifdef ALLOW_USE_MPI
360     IF( mpiMyId .EQ. 0 ) THEN
361     #endif
362     call WRITE_NC_MaskAreaBathy(
363     & 'ECCO','MIT GCM Release 1',
364     & Nx,Ny,Nr,PTRACERS_num,
365     & lon,bounds_lon,lat,bounds_lat,depth,bounds_depth,
366     & MASK,AREA,BATHY,REGION_MASK,REGION_AREA)
367     #ifdef ALLOW_USE_MPI
368     ENDIF
369     #endif
370     _END_MASTER( myThid )
371    
372 dimitri 1.5 #endif /* OCEAN_INVERSION_NETCDF */
373 dimitri 1.6
374 dimitri 1.1 #endif /* ALLOW_PTRACERS */
375    
376     RETURN
377     END

  ViewVC Help
Powered by ViewVC 1.1.22