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/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/TOP/C14 – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/TOP/C14/trcatm_c14.F90 @ 13553

Last change on this file since 13553 was 13553, checked in by hadcv, 3 years ago

Merge in trunk up to [13550]

  • 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( 1, 1, 1, 1 )                 ! from C14b package
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   !!======================================================================
302END MODULE trcatm_c14
Note: See TracBrowser for help on using the repository browser.