Context Navigation

trcatm_c14.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

File size: 13.3 KB

Line
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 + ( rdt / ( 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

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/C14/trcatm_c14.F90 @ 12340

Download in other formats: