1 | MODULE trcatm_c14 |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE trcatm_c14 *** |
---|
4 | !! TOP : read and manages atmospheric values for radiocarbon model |
---|
5 | !!===================================================================== |
---|
6 | !! History: Based on trcini_c14b & trcsms_c14b : |
---|
7 | !! Anne Mouchet |
---|
8 | !!---------------------------------------------------------------------- |
---|
9 | !! trc_atm_c14_ini : initialize c14atm & pco2atm |
---|
10 | !! trc_atm_c14 : read and time interpolate c14atm & pco2atm |
---|
11 | !!---------------------------------------------------------------------- |
---|
12 | USE par_trc ! passive tracers parameters |
---|
13 | USE oce_trc ! shared variables between ocean and passive tracers |
---|
14 | USE trc ! passive tracers common variables |
---|
15 | USE sms_c14 ! c14 simulation type, atm default values... |
---|
16 | |
---|
17 | IMPLICIT NONE |
---|
18 | PRIVATE |
---|
19 | |
---|
20 | PUBLIC trc_atm_c14 ! called in trcsms_c14.F90 |
---|
21 | PUBLIC trc_atm_c14_ini ! called in trcini_c14.F90 |
---|
22 | ! |
---|
23 | !! * Substitutions |
---|
24 | # include "do_loop_substitute.h90" |
---|
25 | !!---------------------------------------------------------------------- |
---|
26 | !! NEMO/TOP 4.0 , NEMO Consortium (2018) |
---|
27 | !! $Id$ |
---|
28 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
29 | !!---------------------------------------------------------------------- |
---|
30 | CONTAINS |
---|
31 | |
---|
32 | SUBROUTINE trc_atm_c14_ini |
---|
33 | !!---------------------------------------------------------------------- |
---|
34 | !! *** ROUTINE trc_c14_ini *** |
---|
35 | !! |
---|
36 | !! ** Purpose : initialisation of sbc for radiocarbon |
---|
37 | !! |
---|
38 | !! ** Method : |
---|
39 | !!---------------------------------------------------------------------- |
---|
40 | ! |
---|
41 | CHARACTER (len=20) :: clfile ! forcing file name |
---|
42 | INTEGER :: ji,jj,jn ! dummy loop indice |
---|
43 | INTEGER :: ierr1,ierr2,ierr3,izco2 ! temporary integers |
---|
44 | INTEGER :: inum1,inum2,incom,iyear ! temporary integers |
---|
45 | REAL(wp) :: ys40 = -40. ! 40 degrees south |
---|
46 | REAL(wp) :: ys20 = -20. ! 20 degrees south |
---|
47 | REAL(wp) :: yn20 = 20. ! 20 degrees north |
---|
48 | REAL(wp) :: yn40 = 40. ! 40 degrees north |
---|
49 | REAL(wp), ALLOCATABLE, DIMENSION(:) :: zco2, zyrco2 ! temporary arrays for swap |
---|
50 | ! |
---|
51 | !!---------------------------------------------------------------------- |
---|
52 | ! |
---|
53 | IF( ln_timing ) CALL timing_start('trc_atm_c14_ini') |
---|
54 | ! |
---|
55 | |
---|
56 | IF( lwp ) WRITE(numout,*) ' ' |
---|
57 | IF( lwp ) WRITE(numout,*) ' trc_atm_c14_ini : initialize atm CO2 & C14-ratio ' |
---|
58 | IF( lwp ) WRITE(numout,*) ' ' |
---|
59 | ! |
---|
60 | tyrc14_now = 0._wp ! initialize |
---|
61 | ! |
---|
62 | IF(kc14typ >= 1) THEN ! Transient atmospheric forcing: CO2 |
---|
63 | ! |
---|
64 | clfile = TRIM( cfileco2 ) |
---|
65 | IF(lwp) WRITE(numout,*) 'Read CO2 atmospheric concentrations file ',clfile |
---|
66 | CALL ctl_opn( inum1, clfile, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) |
---|
67 | REWIND(inum1) |
---|
68 | |
---|
69 | READ(inum1,*) nrecco2,incom |
---|
70 | DO jn = 1, incom ! Skip over descriptor lines |
---|
71 | READ(inum1,'(1x)') |
---|
72 | END DO |
---|
73 | ALLOCATE( spco2(nrecco2), tyrco2(nrecco2) , STAT=ierr1 ) |
---|
74 | IF( ierr1 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate co2 arrays' ) |
---|
75 | ! get CO2 data |
---|
76 | DO jn = 1, nrecco2 |
---|
77 | READ(inum1, *) tyrco2(jn), spco2(jn) |
---|
78 | END DO |
---|
79 | CLOSE(inum1) |
---|
80 | ! |
---|
81 | IF(kc14typ==2) THEN |
---|
82 | ALLOCATE( zco2(nrecco2), zyrco2(nrecco2) ) |
---|
83 | zco2(:)=spco2(:) |
---|
84 | zyrco2(:)=tyrco2(:) |
---|
85 | ! Set CO2 times on AD time scale & swap records : In CO2 file : youngest first |
---|
86 | DO jn = 1, nrecco2 |
---|
87 | izco2=nrecco2-jn+1 |
---|
88 | spco2(izco2)=zco2(jn) |
---|
89 | tyrco2(izco2)=1950._wp-zyrco2(jn) ! BP to AD dates |
---|
90 | END DO |
---|
91 | DEALLOCATE( zco2,zyrco2 ) |
---|
92 | ENDIF |
---|
93 | ! |
---|
94 | ! ! Transient atmospheric forcing: Bomb C14 & Paleo C14 : open file |
---|
95 | ! |
---|
96 | clfile = TRIM( cfilec14 ) |
---|
97 | IF (lwp) WRITE(numout,*) 'Read C-14 atmospheric concentrations file ',clfile |
---|
98 | CALL ctl_opn( inum2, clfile, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) |
---|
99 | REWIND(inum2) |
---|
100 | ! |
---|
101 | ! Bomb C14: 3 zones for atm C14 ! |
---|
102 | IF(kc14typ == 1) THEN ! Transient atmospheric forcing: Bomb C14 |
---|
103 | ! |
---|
104 | READ(inum2,*) nrecc14,incom |
---|
105 | DO jn = 1, incom ! Skip over descriptor lines |
---|
106 | READ(inum2,'(1x)') |
---|
107 | END DO |
---|
108 | ALLOCATE( bomb(nrecc14,nc14zon), tyrc14(nrecc14) , STAT=ierr2 ) |
---|
109 | IF( ierr2 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate c14 arrays' ) |
---|
110 | ! get bomb c14 data |
---|
111 | DO jn = 1, nrecc14 |
---|
112 | READ(inum2,*) tyrc14(jn), bomb(jn,1), bomb(jn,2), bomb(jn,3) |
---|
113 | END DO |
---|
114 | CLOSE(inum2) |
---|
115 | |
---|
116 | ! Linear interpolation of the C-14 source fonction |
---|
117 | ! in linear latitude bands (20N,40N) and (20S,40S) |
---|
118 | !------------------------------------------------------ |
---|
119 | ALLOCATE( fareaz (jpi,jpj ,nc14zon) , STAT=ierr3 ) |
---|
120 | IF( ierr3 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' ) |
---|
121 | ! |
---|
122 | DO_2D_11_11 |
---|
123 | IF( gphit(ji,jj) >= yn40 ) THEN |
---|
124 | fareaz(ji,jj,1) = 0. |
---|
125 | fareaz(ji,jj,2) = 0. |
---|
126 | fareaz(ji,jj,3) = 1. |
---|
127 | ELSE IF( gphit(ji,jj ) <= ys40) THEN |
---|
128 | fareaz(ji,jj,1) = 1. |
---|
129 | fareaz(ji,jj,2) = 0. |
---|
130 | fareaz(ji,jj,3) = 0. |
---|
131 | ELSE IF( gphit(ji,jj) >= yn20 ) THEN |
---|
132 | fareaz(ji,jj,1) = 0. |
---|
133 | fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / yn40 ) |
---|
134 | fareaz(ji,jj,3) = 2. * gphit(ji,jj) / yn40 - 1. |
---|
135 | ELSE IF( gphit(ji,jj) <= ys20 ) THEN |
---|
136 | fareaz(ji,jj,1) = 2. * gphit(ji,jj) / ys40 - 1. |
---|
137 | fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / ys40 ) |
---|
138 | fareaz(ji,jj,3) = 0. |
---|
139 | ELSE |
---|
140 | fareaz(ji,jj,1) = 0. |
---|
141 | fareaz(ji,jj,2) = 1. |
---|
142 | fareaz(ji,jj,3) = 0. |
---|
143 | ENDIF |
---|
144 | END_2D |
---|
145 | ! |
---|
146 | ENDIF |
---|
147 | ! |
---|
148 | ! Paleo C14: 1 zone for atm C14 ! |
---|
149 | IF(kc14typ == 2) THEN ! Transient atmospheric forcing: Paleo C14 |
---|
150 | ! |
---|
151 | READ(inum2,*) nrecc14,incom |
---|
152 | DO jn = 1, incom ! Skip over descriptor lines |
---|
153 | READ(inum2,'(1x)') |
---|
154 | END DO |
---|
155 | ALLOCATE( atmc14(nrecc14), tyrc14(nrecc14) , STAT=ierr2 ) |
---|
156 | IF( ierr2 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate c14 arrays' ) |
---|
157 | ! get past c14 data |
---|
158 | DO jn = 1, nrecc14 |
---|
159 | READ(inum2,*) iyear,incom,incom,atmc14(jn) |
---|
160 | tyrc14(jn)=1950._wp-REAL(iyear,wp) ! BP to AD dates |
---|
161 | END DO |
---|
162 | CLOSE(inum2) |
---|
163 | ! |
---|
164 | ENDIF |
---|
165 | ! |
---|
166 | ! Note on dates: |
---|
167 | ! In files dates have dimension yr; either AD or BP; if BP date is changed into AD here |
---|
168 | ! When dealing with dates previous to 0. AD one needs to set tyrc14_beg to the actual starting year |
---|
169 | ! Do not forget to appropriately set nn_date0 and nn_rstctl in namelist |
---|
170 | ! AND nn_rsttr in namelist_top if offline run |
---|
171 | ! All details are given in NEMO-C14.pdf report |
---|
172 | ! |
---|
173 | tyrc14_now=nyear ! actual initial yr - Bomb |
---|
174 | if(kc14typ == 2) tyrc14_now=nyear+tyrc14_beg-1 ! actual initial yr - Paleo |
---|
175 | ! ! we suppose we start on tyrc14_now/01/01 @ 0h |
---|
176 | m1_c14= 1 |
---|
177 | m1_co2= 1 |
---|
178 | DO jn = 1,nrecco2 |
---|
179 | IF ( tyrc14_now >= tyrco2(jn) ) m1_co2 = jn ! index of first co2 record to consider |
---|
180 | END DO |
---|
181 | DO jn = 1,nrecc14 |
---|
182 | IF ( tyrc14_now >= tyrc14(jn) ) m1_c14 = jn ! index of first c14 record to consider |
---|
183 | END DO |
---|
184 | IF (lwp) WRITE(numout,*) 'Initial yr for experiment', tyrc14_now |
---|
185 | IF (lwp) WRITE(numout,*) ' CO2 & C14 start years:', tyrco2(m1_co2),tyrc14(m1_c14) |
---|
186 | ! |
---|
187 | m2_c14= m1_c14 |
---|
188 | m2_co2= m1_co2 |
---|
189 | ! |
---|
190 | ENDIF |
---|
191 | ! |
---|
192 | IF( ln_timing ) CALL timing_stop('trc_atm_c14_ini') |
---|
193 | ! |
---|
194 | END SUBROUTINE trc_atm_c14_ini |
---|
195 | |
---|
196 | |
---|
197 | SUBROUTINE trc_atm_c14( kt, co2sbc, c14sbc ) |
---|
198 | !!---------------------------------------------------------------------- |
---|
199 | !! *** ROUTINE trc_flx *** |
---|
200 | !! |
---|
201 | !! ** Purpose : provides sbc for co2 & c14 at kt |
---|
202 | !! |
---|
203 | !! ** Method : read files |
---|
204 | !! |
---|
205 | !! ** Action : atmospheric values interpolated at time-step kt |
---|
206 | !!---------------------------------------------------------------------- |
---|
207 | INTEGER , INTENT(in ) :: kt ! ocean time-step |
---|
208 | REAL(wp), DIMENSION(:,:), INTENT( out) :: c14sbc ! atm c14 ratio |
---|
209 | REAL(wp), INTENT( out) :: co2sbc ! atm co2 p |
---|
210 | INTEGER :: jz ! dummy loop indice |
---|
211 | REAL(wp) :: zdint,zint ! work |
---|
212 | REAL(wp), DIMENSION(nc14zon) :: zonbc14 ! work |
---|
213 | ! |
---|
214 | !!---------------------------------------------------------------------- |
---|
215 | ! |
---|
216 | IF( ln_timing ) CALL timing_start('trc_atm_c14') |
---|
217 | ! |
---|
218 | IF( kc14typ == 0) THEN |
---|
219 | co2sbc=pco2at |
---|
220 | c14sbc(:,:)=rc14at |
---|
221 | ENDIF |
---|
222 | ! |
---|
223 | IF(kc14typ >= 1) THEN ! Transient C14 & CO2 |
---|
224 | ! |
---|
225 | tyrc14_now = tyrc14_now + ( rn_Dt / ( rday * nyear_len(1)) ) ! current time step in yr relative to tyrc14_beg |
---|
226 | ! |
---|
227 | ! CO2 -------------------------------------------------------- |
---|
228 | ! |
---|
229 | ! time interpolation of CO2 concentrations ! if out of record keeps first/last value |
---|
230 | IF( tyrc14_now > tyrco2(m2_co2) ) THEN ! next interval |
---|
231 | m1_co2 = m2_co2 |
---|
232 | m2_co2 = MIN ( m2_co2 + 1 , nrecco2 ) |
---|
233 | ENDIF |
---|
234 | ! |
---|
235 | zdint = tyrco2(m2_co2) - tyrco2(m1_co2) |
---|
236 | co2sbc = spco2(m2_co2) ! if out of record keeps first/last value |
---|
237 | zint = 0._wp |
---|
238 | IF ( zdint > 0._wp ) THEN ! if within record interpolate: |
---|
239 | zint = ( tyrco2(m2_co2) - tyrc14_now ) / zdint |
---|
240 | co2sbc = spco2(m2_co2) + zint * ( spco2(m1_co2) - spco2(m2_co2) ) |
---|
241 | ENDIF |
---|
242 | ! |
---|
243 | IF( lwp .AND. kt == nitend ) THEN |
---|
244 | WRITE(numout, '(3(A,F12.4))') 't1/tn/t2:',tyrco2(m1_co2),'/', tyrc14_now,'/',tyrco2(m2_co2) |
---|
245 | WRITE(numout, *) 'CO2:',spco2(m1_co2),co2sbc ,spco2(m2_co2) |
---|
246 | ENDIF |
---|
247 | ! |
---|
248 | ! C14 -------------------------------------------------------- |
---|
249 | ! |
---|
250 | ! time interpolation of C14 concentrations |
---|
251 | IF ( tyrc14_now > tyrc14(m2_c14) ) THEN ! next interval |
---|
252 | m1_c14 = m2_c14 |
---|
253 | m2_c14 = MIN ( m2_c14 + 1 , nrecc14 ) |
---|
254 | ENDIF |
---|
255 | zdint = tyrc14(m2_c14) - tyrc14(m1_c14) |
---|
256 | zint=0._wp |
---|
257 | IF ( zdint > 0._wp ) zint = ( tyrc14(m2_c14) - tyrc14_now ) / zdint ! if within record |
---|
258 | IF( lwp .AND. kt == nitend ) & |
---|
259 | & WRITE(numout,'(3(A,F12.4))') 't1/tn/t2:',tyrc14(m1_c14),'/', tyrc14_now,'/',tyrc14(m2_c14) |
---|
260 | ! |
---|
261 | ! ------- Bomb C14 ------------------------------------------ |
---|
262 | ! |
---|
263 | IF( kc14typ == 1) THEN |
---|
264 | ! ! time interpolation |
---|
265 | zonbc14(:) = bomb(m2_c14,:) ! if out of record keeps first/last value |
---|
266 | ! ! if within record interpolate: |
---|
267 | IF ( zdint > 0._wp ) zonbc14(:) = bomb(m2_c14,:) + zint * ( bomb(m1_c14,:) - bomb(m2_c14,:) ) |
---|
268 | ! |
---|
269 | IF(lwp .AND. kt == nitend ) & |
---|
270 | & WRITE(numout, *) 'C14:',bomb(m1_c14,1),zonbc14(1),bomb(m2_c14,1) |
---|
271 | ! Transform DeltaC14 --> C14 ratio |
---|
272 | zonbc14(:) = 1._wp + zonbc14(:)/1.d03 |
---|
273 | ! |
---|
274 | ! For each (i,j)-box, with information from the fractional area |
---|
275 | ! (zonmean), computes area-weighted mean to give the atmospheric C-14 |
---|
276 | ! ---------------------------------------------------------------- |
---|
277 | c14sbc(:,:) = zonbc14(1) * fareaz(:,:,1) & |
---|
278 | & + zonbc14(2) * fareaz(:,:,2) & |
---|
279 | & + zonbc14(3) * fareaz(:,:,3) |
---|
280 | ENDIF |
---|
281 | ! |
---|
282 | ! ------- Paleo C14 ----------------------------------------- |
---|
283 | ! |
---|
284 | IF( kc14typ == 2 ) THEN |
---|
285 | ! ! time interpolation |
---|
286 | zonbc14(1) = atmc14(m2_c14) ! if out of record keeps first/last value |
---|
287 | ! ! if within record interpolate: |
---|
288 | IF ( zdint > 0._wp ) zonbc14(1) = atmc14(m2_c14) + zint * ( atmc14(m1_c14) - atmc14(m2_c14) ) |
---|
289 | IF(lwp .AND. kt == nitend ) & |
---|
290 | & WRITE(numout, *) 'C14: ',atmc14(m1_c14),zonbc14(1),atmc14(m2_c14) |
---|
291 | ! Transform DeltaC14 --> C14 ratio |
---|
292 | c14sbc(:,:) = 1._wp + zonbc14(1)/1.d03 |
---|
293 | ENDIF |
---|
294 | ! |
---|
295 | ENDIF |
---|
296 | ! |
---|
297 | IF( ln_timing ) CALL timing_stop('trc_atm_c14') |
---|
298 | ! |
---|
299 | END SUBROUTINE trc_atm_c14 |
---|
300 | |
---|
301 | !!====================================================================== |
---|
302 | END MODULE trcatm_c14 |
---|