/[MITgcm]/MITgcm_contrib/ocean_inversion_project/code_bombC14/ini_tr1.F
ViewVC logotype

Contents of /MITgcm_contrib/ocean_inversion_project/code_bombC14/ini_tr1.F

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


Revision 1.3 - (show annotations) (download)
Thu May 25 20:57:50 2006 UTC (19 years, 2 months ago) by dimitri
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +47 -2 lines
preparing bomb C14 experiment

1 C $Header: /u/gcmpack/MITgcm_contrib/ocean_inversion_project/code_bombC14/ini_tr1.F,v 1.2 2006/05/25 17:16:56 dimitri Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: INI_TR1
8 C !INTERFACE:
9 SUBROUTINE INI_TR1( myThid )
10 C !DESCRIPTION: \bv
11 C *==========================================================*
12 C | SUBROUTINE INI_TR1
13 C | o Set initial tracer 1 distribution.
14 C *==========================================================*
15 C | Passive tracers 1-N can be initialised so as to dye fluid.
16 C | This routine is hardcoded for Nir's bomb C14 experiment.
17 C *==========================================================*
18 C \ev
19
20 C !USES:
21 IMPLICIT NONE
22 C === Global variables ===
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "GRID.h"
27 #include "DYNVARS.h"
28 #ifdef ALLOW_PASSIVE_TRACER
29 #include "TR1.h"
30 #endif
31
32 C !INPUT/OUTPUT PARAMETERS:
33 C == Routine arguments ==
34 C myThid - Number of this instance of INI_TR1
35 INTEGER myThid
36
37 #ifdef ALLOW_PASSIVE_TRACER
38
39 C !LOCAL VARIABLES:
40 C == Local variables ==
41 C iC, jC - Center of domain
42 C iD, jD - Disitance from domain center.
43 C rad - Radius of initial patch
44 C rD - Radial displacement of point I,J
45 C iG, jG - Global coordinate index
46 C bi,bj - Loop counters
47 C I,J,K
48 INTEGER iC, jC, iD, jD
49 INTEGER iG, jG
50 INTEGER bi, bj
51 INTEGER I, J, K, localWarnings
52 _RL rad, rD
53 CHARACTER*(MAX_LEN_MBUF) msgBuf
54 CEOP
55
56 _BARRIER
57
58 C-- Initialise tracer to inline.
59 DO bj = myByLo(myThid), myByHi(myThid)
60 DO bi = myBxLo(myThid), myBxHi(myThid)
61 DO K=1,Nr
62 DO J=1,sNy
63 DO I=1,sNx
64 tr1(I,J,K,bi,bj) = 0. _d 0
65 ENDDO
66 ENDDO
67 ENDDO
68 ENDDO
69 ENDDO
70 C Set initial tendency terms
71 localWarnings=0
72 DO bj = myByLo(myThid), myByHi(myThid)
73 DO bi = myBxLo(myThid), myBxHi(myThid)
74 DO K=1,Nr
75 DO J=1,sNy
76 DO I=1,sNx
77 gtr1 (I,J,K,bi,bj) = 0. _d 0
78 gtr1NM1(I,J,K,bi,bj) = 0. _d 0
79 ENDDO
80 ENDDO
81 ENDDO
82 ENDDO
83 ENDDO
84 C
85 _EXCH_XYZ_R8(tr1 , myThid )
86 _EXCH_XYZ_R8(gtr1 , myThid )
87 _EXCH_XYZ_R8(gtr1NM1 , myThid )
88
89 C-- Input file names
90 FiceFile = 'FICE'
91 XkwFile = 'XKW'
92 PatmFile = 'PATM'
93 WOASssFile = 'WOA01_S'
94 WOASstFile = 'WOA01_T'
95 SurfDicPreindFile = 'surf_dic_preind'
96 SurfDicAnthroFile = 'surf_dic_anthro'
97 SurfNatD14CFile = 'surface_naturalD14C'
98
99 C-- Compute numbers needed to determine DC14 and CO2_atm
100 CALL read_dc14
101 CALL read_CO2
102
103 _BEGIN_MASTER(myThid)
104 #define OCMIP_USE_INTERP
105 C-- OCMIP_USE_INTERP option is useful for high-resolution integrations.
106 C For low-resolution, less that 1-deg, it's best to generate files
107 C separately because of sea-ice fraction.
108 C Caution: OCMIP_USE_INTERP as used is not thread-safe.
109 #ifdef OCMIP_USE_INTERP
110 CALL EXF_INTERP( SurfDicPreindFile,readBinaryPrec,
111 & SurfDicPreind,1,
112 & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
113 CALL EXF_INTERP( SurfDicAnthroFile,readBinaryPrec,
114 & SurfDicAnthro,1,
115 & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
116 CALL EXF_INTERP( SurfNatD14CFile,readBinaryPrec,
117 & SurfNatD14C,1,
118 & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
119
120 #else
121 CALL MDSREADFIELD ( SurfDicPreindFile, readBinaryPrec,
122 & 'RS', 1, SurfDicPreind, 1,, myThid )
123 CALL MDSREADFIELD ( SurfDicAnthroFile, readBinaryPrec,
124 & 'RS', 1, SurfDicAnthro, 1,, myThid )
125 CALL MDSREADFIELD ( SurfNatD14CFile, readBinaryPrec,
126 & 'RS', 1, SurfNatD14C, 1,, myThid )
127 #endif
128 _END_MASTER(myThid)
129 _EXCH_XY_R4(SurfDicPreind , myThid )
130 _EXCH_XY_R4(SurfDicAnthro , myThid )
131 _EXCH_XY_R4(SurfNatD14C , myThid )
132
133 #endif /* ALLOW_PASSIVE_TRACER */
134
135 RETURN
136 END
137
138 C====================================================================
139
140 SUBROUTINE read_dc14
141
142 C --------------------------------------------------------------------
143 C Reads temporal history of atmospheric DC14 in permil
144 C
145 C DC14_atm_year = year of measurement (mid-year)
146 C DC14_atm_n = DC14 in northern hemisphere
147 C DC14_atm_eq = DC14 at equator
148 C DC14_atm_s = DC14 in southern hemisphere
149 C --------------------------------------------------------------------
150
151 IMPLICIT NONE
152 #include "SIZE.h"
153 #include "EEPARAMS.h"
154 #include "TR1.h"
155
156 INTEGER n, iu
157 PARAMETER (iu=10)
158
159 CHARACTER*80 filen
160 C
161 C OPEN FILE
162 C ---------
163 filen='DC14_atm'
164 OPEN(UNIT=iu,FILE=filen,STATUS='old')
165 C
166 C READ FILE
167 C ---------
168 DO n = 1, 240
169 READ(iu,*)
170 + DC14_atm_year(n), DC14_atm_n(n),
171 + DC14_atm_eq (n), DC14_atm_s(n)
172 C WRITE(6,100)
173 C + DC14_atm_year(n), DC14_atm_n(n),
174 C + DC14_atm_eq (n), DC14_atm_s(n)
175 END DO
176
177 100 FORMAT(f7.2, 4f8.2)
178
179 RETURN
180 END
181
182 C====================================================================
183
184 SUBROUTINE read_CO2
185
186 C --------------------------------------------------------------------
187 C Reads temporal history of atmospheric CO2 in ppm
188 C
189 C CO2_atm_year = year of measurement
190 C CO2_atm = atmospheric CO2
191 C --------------------------------------------------------------------
192
193 IMPLICIT NONE
194 #include "SIZE.h"
195 #include "EEPARAMS.h"
196 #include "TR1.h"
197
198 INTEGER n, iu
199 PARAMETER (iu=10)
200
201 CHARACTER*80 filen
202 C
203 C OPEN FILE
204 C ---------
205 filen='CO2_atm'
206 OPEN(UNIT=iu,FILE=filen,STATUS='old')
207 C
208 C READ FILE
209 C ---------
210 DO n = 1, 478
211 READ(iu,*) CO2_atm_year(n), CO2_atm(n)
212 C WRITE(6,100) CO2_atm_year(n), CO2_atm(n)
213 END DO
214
215 100 FORMAT(f7.2, f8.2)
216
217 RETURN
218 END

  ViewVC Help
Powered by ViewVC 1.1.22