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

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

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

revision 1.1 by dimitri, Thu Sep 18 02:33:38 2003 UTC revision 1.9 by dimitri, Wed Oct 6 20:26:54 2004 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "PTRACERS_OPTIONS.h"  #include "PTRACERS_OPTIONS.h"
 #ifdef ALLOW_GCHEM  
 # include "GCHEM_OPTIONS.h"  
 #endif  
5    
6  CBOP  CBOP
7  C !ROUTINE: PTRACERS_INIT  C !ROUTINE: PTRACERS_INIT
# Line 14  C !INTERFACE: ========================== Line 11  C !INTERFACE: ==========================
11    
12  C !DESCRIPTION:  C !DESCRIPTION:
13  C     Initialize PTRACERS data structures  C     Initialize PTRACERS data structures
14  cdm   This file is customized to compute CO2 perturbations from 30 ocean  C     This file is customized to compute CO2 perturbations from
15  cdm   regions for Gruber's ocean inversion project.  C     30 ocean regions for Gruber's ocean inversion project.
16    
17  C !USES: ===============================================================  C !USES: ===============================================================
18        IMPLICIT NONE        IMPLICIT NONE
19  #include "SIZE.h"  #include "SIZE.h"
20  #include "EEPARAMS.h"  #include "EEPARAMS.h"
21    #include "EESUPPORT.h"
22  #include "PARAMS.h"  #include "PARAMS.h"
23  #include "GRID.h"  #include "GRID.h"
24  #include "PTRACERS.h"  #include "PTRACERS.h"
 cswdptr -- add ---  
 #ifdef ALLOW_GCHEM  
 # include "GCHEM.h"  
 #endif  
 cswdptr --- end add --  
25    
26    
27  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
# Line 42  C  none Line 35  C  none
35    
36  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
37  C  i,j,k,bi,bj,iTracer  :: loop indices  C  i,j,k,bi,bj,iTracer  :: loop indices
38        INTEGER i,j,k,bi,bj,iTracer        INTEGER i,j,k,bi,bj,iTracer,iMonth,iUnit,iRec
39        CHARACTER*(10) suff        CHARACTER*(10) suff
40  #ifndef ALLOW_GCHEM        _RL SumPtracer
41        INTEGER tIter0  
42        PARAMETER ( tIter0 = 0 )  #ifdef OCEAN_INVERSION_NETCDF
43  #endif        Real*8 global(Nx,Ny)
44          _RL    local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45          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  CEOP  CEOP
56                    
57  C Loop over tracers  C Loop over tracers
58        DO iTracer = 1, PTRACERS_num        DO iTracer = 1, PTRACERS_num
59    
60  C Loop over tiles  C Loop over tiles
61        DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
62         DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
63    
64  C Initialize arrays in common blocks :  C Initialize arrays in common blocks :
65          DO k=1,Nr           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              ENDDO
73             ENDDO
74           DO j=1-Oly,sNy+OLy           DO j=1-Oly,sNy+OLy
75            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
76             pTracer(i,j,k,bi,bj,iTracer) = 0. _d 0             surfaceTendencyPtr(i,j,bi,bj,iTracer) = 0. _d 0
            gPtr(i,j,k,bi,bj,iTracer)    = 0. _d 0  
            gPtrNM1(i,j,k,bi,bj,iTracer) = 0. _d 0  
77            ENDDO            ENDDO
78           ENDDO           ENDDO
         ENDDO  
79    
80  C end bi,bj loops  C end bi,bj loops
81            ENDDO
82         ENDDO         ENDDO
       ENDDO  
83    
84  C end of Tracer loop  C end of Tracer loop
85        ENDDO        ENDDO
86    
87  C Read from a pickup file if nIter0  C Read from a pickup file if nIter0 is not zero
88  cswdptr IF (nIter0.NE.0) THEN        IF (nIter0.NE.0) THEN
 cswdptr -- change --  
       IF (nIter0.GT.tIter0) THEN  
89  C--     Suffix for pickup files  C--     Suffix for pickup files
90         IF (pickupSuff.EQ.' ') THEN         IF (pickupSuff.EQ.' ') THEN
91           WRITE(suff,'(I10.10)') nIter0           WRITE(suff,'(I10.10)') nIter0
# Line 88  C--     Suffix for pickup files Line 95  C--     Suffix for pickup files
95         CALL PTRACERS_READ_CHECKPOINT( nIter0,myThid )         CALL PTRACERS_READ_CHECKPOINT( nIter0,myThid )
96        ENDIF        ENDIF
97    
98  cdm   Initialize pTracerMasks  C     Initialize pTracerMasks
99        CALL PTRACERS_READ_MASK ( mythid )        CALL PTRACERS_READ_MASK ( mythid )
100    
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    C     It is assumed that each year is 365.24167 days (31556880 s)
107    C     long and that each month is 2629740 s.
108    
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                SumPtracer = SumPtracer + 2629740. * rA(I,J,bi,bj) *
117         &             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    
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    
151    #ifdef OCEAN_INVERSION_NETCDF
152    
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          DO j=1,Ny
281           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             ENDDO
331            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          ENDDO
340    
341    C--   depth, bounds_depth
342          DO k=1,Nr
343           depth        (k  ) = abs ( rC (k)   )
344           bounds_depth (k,1) = abs ( rF (k+1) )
345           bounds_depth (k,2) = abs ( rF (k)   )
346          ENDDO
347    
348          DO iTracer = 1, PTRACERS_num
349           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          ENDDO
357    
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    #endif /* OCEAN_INVERSION_NETCDF */
374    
375  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
376    
377        RETURN        RETURN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22