/[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.8 by dimitri, Thu Dec 18 06:11:14 2003 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
# Line 24  C !USES: =============================== Line 21  C !USES: ===============================
21  #include "PARAMS.h"  #include "PARAMS.h"
22  #include "GRID.h"  #include "GRID.h"
23  #include "PTRACERS.h"  #include "PTRACERS.h"
 cswdptr -- add ---  
 #ifdef ALLOW_GCHEM  
 # include "GCHEM.h"  
 #endif  
 cswdptr --- end add --  
24    
25    
26  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
# Line 42  C  none Line 34  C  none
34    
35  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
36  C  i,j,k,bi,bj,iTracer  :: loop indices  C  i,j,k,bi,bj,iTracer  :: loop indices
37        INTEGER i,j,k,bi,bj,iTracer        INTEGER i,j,k,bi,bj,iTracer,iMonth,iUnit,iRec
38        CHARACTER*(10) suff        CHARACTER*(10) suff
39  #ifndef ALLOW_GCHEM        _RL SumPtracer
40        INTEGER tIter0  
41        PARAMETER ( tIter0 = 0 )  #ifdef OCEAN_INVERSION_NETCDF
42  #endif        Real*8 global(Nx,Ny)
43          _RL    local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44          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  CEOP  CEOP
55                    
56  C Loop over tracers  C Loop over tracers
57        DO iTracer = 1, PTRACERS_num        DO iTracer = 1, PTRACERS_num
58    
59  C Loop over tiles  C Loop over tiles
60        DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
61         DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
62    
63  C Initialize arrays in common blocks :  C Initialize arrays in common blocks :
64          DO k=1,Nr           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              ENDDO
72             ENDDO
73           DO j=1-Oly,sNy+OLy           DO j=1-Oly,sNy+OLy
74            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
75             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  
76            ENDDO            ENDDO
77           ENDDO           ENDDO
         ENDDO  
78    
79  C end bi,bj loops  C end bi,bj loops
80            ENDDO
81         ENDDO         ENDDO
       ENDDO  
82    
83  C end of Tracer loop  C end of Tracer loop
84        ENDDO        ENDDO
85    
86  C Read from a pickup file if nIter0  C Read from a pickup file if nIter0 is not zero
87  cswdptr IF (nIter0.NE.0) THEN        IF (nIter0.NE.0) THEN
 cswdptr -- change --  
       IF (nIter0.GT.tIter0) THEN  
88  C--     Suffix for pickup files  C--     Suffix for pickup files
89         IF (pickupSuff.EQ.' ') THEN         IF (pickupSuff.EQ.' ') THEN
90           WRITE(suff,'(I10.10)') nIter0           WRITE(suff,'(I10.10)') nIter0
# Line 88  C--     Suffix for pickup files Line 94  C--     Suffix for pickup files
94         CALL PTRACERS_READ_CHECKPOINT( nIter0,myThid )         CALL PTRACERS_READ_CHECKPOINT( nIter0,myThid )
95        ENDIF        ENDIF
96    
97  cdm   Initialize pTracerMasks  C     Initialize pTracerMasks
98        CALL PTRACERS_READ_MASK ( mythid )        CALL PTRACERS_READ_MASK ( mythid )
99    
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.24167 days (31556880 s)
106    C     long and that each month is 2629740 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 + 2629740. * 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    
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    
150    #ifdef OCEAN_INVERSION_NETCDF
151    
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          DO j=1,Ny
280           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             ENDDO
330            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          ENDDO
339    
340    C--   depth, bounds_depth
341          DO k=1,Nr
342           depth        (k  ) = abs ( rC (k)   )
343           bounds_depth (k,1) = abs ( rF (k+1) )
344           bounds_depth (k,2) = abs ( rF (k)   )
345          ENDDO
346    
347          DO iTracer = 1, PTRACERS_num
348           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          ENDDO
356    
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    #endif /* OCEAN_INVERSION_NETCDF */
373    
374  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
375    
376        RETURN        RETURN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22