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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trcatm_c14.F90 in NEMO/trunk/src/TOP/C14 – NEMO

source: NEMO/trunk/src/TOP/C14/trcatm_c14.F90 @ 12377

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

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 13.3 KB
Line 
1MODULE 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   !!----------------------------------------------------------------------
30CONTAINS
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   !!======================================================================
302END MODULE trcatm_c14
Note: See TracBrowser for help on using the repository browser.