source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/C14/trcatm_c14.F90 @ 10372

Last change on this file since 10372 was 10069, checked in by nicolasmartin, 2 years ago

Fix mistakes of previous commit on SVN keywords property

  • Property svn:keywords set to Id
File size: 13.4 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   !!----------------------------------------------------------------------
24   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
25   !! $Id$
26   !! Software governed by the CeCILL license (see ./LICENSE)
27   !!----------------------------------------------------------------------
28CONTAINS
29
30   SUBROUTINE trc_atm_c14_ini
31      !!----------------------------------------------------------------------
32      !!                   ***  ROUTINE trc_c14_ini  ***
33      !!                   
34      !! ** Purpose :   initialisation of sbc for radiocarbon
35      !!
36      !! ** Method  :
37      !!----------------------------------------------------------------------
38      !
39      CHARACTER (len=20)        :: clfile                        ! forcing file name
40      INTEGER                   :: ji,jj,jn                        ! dummy loop indice
41      INTEGER                   :: ierr1,ierr2,ierr3,izco2         ! temporary integers
42      INTEGER                   :: inum1,inum2,incom,iyear         ! temporary integers
43      REAL(wp) ::   ys40 = -40.     ! 40 degrees south
44      REAL(wp) ::   ys20 = -20.     ! 20 degrees south
45      REAL(wp) ::   yn20 =  20.     ! 20 degrees north
46      REAL(wp) ::   yn40 =  40.     ! 40 degrees north
47      REAL(wp), ALLOCATABLE, DIMENSION(:)   :: zco2, zyrco2  ! temporary arrays for swap
48      !
49      !!----------------------------------------------------------------------
50      !
51      IF( ln_timing )   CALL timing_start('trc_atm_c14_ini')
52      !
53     
54      IF( lwp ) WRITE(numout,*) ' '
55      IF( lwp ) WRITE(numout,*) ' trc_atm_c14_ini : initialize atm CO2 & C14-ratio '
56      IF( lwp ) WRITE(numout,*) ' '
57      !
58      tyrc14_now = 0._wp   ! initialize
59      !
60      IF(kc14typ >= 1) THEN  ! Transient atmospheric forcing: CO2
61      !
62         clfile = TRIM( cfileco2 )
63         IF(lwp) WRITE(numout,*) 'Read CO2 atmospheric concentrations file ',clfile
64         CALL ctl_opn( inum1, clfile, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
65         REWIND(inum1)
66     
67         READ(inum1,*) nrecco2,incom
68         DO jn = 1, incom        ! Skip over descriptor lines
69            READ(inum1,'(1x)')
70         END DO
71         ALLOCATE( spco2(nrecco2), tyrco2(nrecco2)       , STAT=ierr1 )
72         IF( ierr1 /= 0 )   CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate co2 arrays' )
73      ! get  CO2 data
74         DO jn = 1, nrecco2
75            READ(inum1, *)  tyrco2(jn), spco2(jn)
76         END DO
77         CLOSE(inum1)
78      !
79         IF(kc14typ==2) THEN
80            ALLOCATE( zco2(nrecco2), zyrco2(nrecco2) )
81            zco2(:)=spco2(:)
82            zyrco2(:)=tyrco2(:)
83      ! Set CO2 times on AD time scale & swap records : In CO2 file : youngest first
84            DO jn = 1, nrecco2
85               izco2=nrecco2-jn+1
86               spco2(izco2)=zco2(jn)
87               tyrco2(izco2)=1950._wp-zyrco2(jn)         ! BP to AD dates
88            END DO
89            DEALLOCATE( zco2,zyrco2 )
90         ENDIF
91      !
92      !        ! Transient atmospheric forcing: Bomb C14 & Paleo C14 : open file
93      !
94         clfile = TRIM( cfilec14 )
95         IF (lwp) WRITE(numout,*) 'Read C-14 atmospheric concentrations file ',clfile
96         CALL ctl_opn( inum2, clfile, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
97         REWIND(inum2)
98      !
99      ! Bomb C14: 3 zones for atm C14 !
100         IF(kc14typ == 1) THEN  ! Transient atmospheric forcing: Bomb C14
101      !
102            READ(inum2,*) nrecc14,incom
103            DO jn = 1, incom        ! Skip over descriptor lines
104               READ(inum2,'(1x)')
105            END DO
106            ALLOCATE( bomb(nrecc14,nc14zon), tyrc14(nrecc14)       , STAT=ierr2 )
107            IF( ierr2 /= 0 )   CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate c14 arrays' )
108      ! get bomb c14 data
109            DO jn = 1, nrecc14
110               READ(inum2,*) tyrc14(jn), bomb(jn,1), bomb(jn,2), bomb(jn,3)
111            END DO
112            CLOSE(inum2)
113
114       ! Linear  interpolation of the C-14 source fonction
115       ! in linear latitude bands  (20N,40N) and (20S,40S)
116       !------------------------------------------------------
117            ALLOCATE( fareaz  (jpi,jpj ,nc14zon) , STAT=ierr3 )
118            IF( ierr3 /= 0 )   CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' )
119      !
120            DO jj = 1 , jpj                       ! from C14b package
121              DO ji = 1 , jpi
122                 IF( gphit(ji,jj) >= yn40 ) THEN
123                    fareaz(ji,jj,1) = 0.
124                    fareaz(ji,jj,2) = 0.
125                    fareaz(ji,jj,3) = 1.
126                 ELSE IF( gphit(ji,jj ) <= ys40) THEN
127                    fareaz(ji,jj,1) = 1.
128                    fareaz(ji,jj,2) = 0.
129                    fareaz(ji,jj,3) = 0.
130                 ELSE IF( gphit(ji,jj) >= yn20 ) THEN
131                    fareaz(ji,jj,1) = 0.
132                    fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / yn40 )
133                    fareaz(ji,jj,3) = 2. * gphit(ji,jj) / yn40 - 1.
134                 ELSE IF( gphit(ji,jj) <= ys20 ) THEN
135                    fareaz(ji,jj,1) = 2. * gphit(ji,jj) / ys40 - 1.
136                    fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / ys40 )
137                    fareaz(ji,jj,3) = 0.
138                 ELSE
139                    fareaz(ji,jj,1) = 0.
140                    fareaz(ji,jj,2) = 1.
141                    fareaz(ji,jj,3) = 0.
142                 ENDIF
143              END DO
144           END DO
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.