/[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.9 - (hide annotations) (download)
Wed Oct 6 20:26:54 2004 UTC (20 years, 9 months ago) by dimitri
Branch: MAIN
CVS Tags: HEAD
Changes since 1.8: +2 -1 lines
Modified eesupp/inc/EEPARAMS.h
Added ocean_inversion_project/*

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

  ViewVC Help
Powered by ViewVC 1.1.22