source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File size: 14.7 KB
Line 
1MODULE trcsms_c14b
2   !!======================================================================
3   !!                      ***  MODULE trcsms_c14b  ***
4   !! TOP : Bomb C14 main module
5   !!======================================================================
6   !! History     -   ! 1994-05 ( J. Orr ) original code
7   !!            1.0  ! 2006-02 ( J.M. Molines )  Free form + modularity
8   !!            2.0  ! 2008-12 ( C. Ethe ) reorganisation
9   !!            4.0  ! 2011-02 ( A.R. Porter, STFC Daresbury ) Dynamic memory
10   !!----------------------------------------------------------------------
11#if defined key_c14b
12   !!----------------------------------------------------------------------
13   !!   'key_c14b'                                         Bomb C14 tracer
14   !!----------------------------------------------------------------------
15   !!   trc_sms_c14b :  compute and add C14 suface forcing to C14 trends
16   !!----------------------------------------------------------------------
17   USE oce_trc       ! Ocean variables
18   USE par_trc       ! TOP parameters
19   USE trc           ! TOP variables
20   USE trd_oce
21   USE trdtrc
22   USE iom           ! I/O library
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   trc_sms_c14b       ! called in trcsms.F90
28   PUBLIC   trc_sms_c14b_alloc ! called in trcini_c14b.F90
29
30   INTEGER , PUBLIC, PARAMETER ::   jpmaxrec  = 240           ! temporal parameter
31   INTEGER , PUBLIC, PARAMETER ::   jpmaxrec2 = 2 * jpmaxrec  !
32   INTEGER , PUBLIC, PARAMETER ::   jpzon     = 3             ! number of zones
33
34   INTEGER , PUBLIC    ::   ndate_beg_b      ! initial calendar date (aammjj) for C14
35   INTEGER , PUBLIC    ::   nyear_beg_b      ! initial year for C14
36   INTEGER , PUBLIC    ::   nyear_res_b      ! restoring time constant (year)
37   INTEGER , PUBLIC    ::   nyear_beg        ! initial year (aa)
38
39   REAL(wp), PUBLIC,                    DIMENSION(jpmaxrec,jpzon) ::   bomb       !: C14 atm data (3 zones)
40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)          ::   fareaz     !: Spatial Interpolation Factors
41   REAL(wp), PUBLIC,                    DIMENSION(jpmaxrec2)      ::   spco2      !: Atmospheric CO2
42 
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)            ::   qtr_c14    !: flux at surface
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)            ::   qint_c14   !: cumulative flux
45
46   REAL(wp) ::   xlambda, xdecay, xaccum       ! C14 decay coef. 
47   REAL(wp) ::   xconv1 = 1._wp                ! conversion from to
48   REAL(wp) ::   xconv2 = 0.01_wp / 3600._wp   ! conversion from cm/h to m/s:
49   REAL(wp) ::   xconv3 = 1.e+3_wp             ! conversion from mol/l/atm to mol/m3/atm
50
51   !! * Substitutions
52#  include "top_substitute.h90"
53
54   !!----------------------------------------------------------------------
55   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE trc_sms_c14b( kt )
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE trc_sms_c14b  ***
64      !!
65      !! ** Purpose :   Compute the surface boundary contition on C14bomb
66      !!      passive tracer associated with air-mer fluxes and add it to
67      !!      the general trend of tracers equations.
68      !!
69      !! ** Original comments from J. Orr :
70      !!
71      !!      Calculates the input of Bomb C-14 to the surface layer of OPA
72      !!
73      !!      James Orr, LMCE, 28 October 1992
74      !!
75      !!      Initial approach mimics that of Toggweiler, Dixon, & Bryan (1989)
76      !!      (hereafter referred to as TDB) with constant gas exchange,
77      !!      although in this case, a perturbation approach is used for
78      !!      bomb-C14 so that both the ocean and atmosphere begin at zero.
79      !!      This saves tremendous amounts of computer time since no
80      !!      equilibrum run is first required (i.e., for natural C-14).
81      !!      Note: Many sensitivity tests can be run with this approach and
82      !!            one never has to make a run for natural C-14; otherwise,
83      !!            a run for natural C-14 must be run each time that one
84      !!            changes a model parameter!
85      !!
86      !!
87      !!      19 August 1993: Modified to consider Atmospheric C-14 fom IPCC.
88      !!      That is, the IPCC has provided a C-14 atmospheric record (courtesy
89      !!      of Martin Heimann) for model calibration.  This model spans from
90      !!      preindustrial times to present, in a format different than that
91      !!      given by TDB.  It must be converted to the ORR C-14 units used
92      !!      here, although in this case, the perturbation includes not only
93      !!      bomb C-14 but changes due to the Suess effect.
94      !!
95      !!----------------------------------------------------------------------
96      !
97      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
98      !
99      INTEGER :: ji, jj, jk, jz     ! dummy loop indices
100      INTEGER :: iyear_beg , iyear_beg1, iyear_end1 
101      INTEGER :: iyear_beg2, iyear_end2 
102      INTEGER :: imonth1, im1, in1 
103      INTEGER :: imonth2, im2, in2 
104      REAL(wp), DIMENSION(jpzon) :: zonbc14       !: time interp atm C14
105      REAL(wp)                   :: zpco2at       !: time interp atm C02
106      REAL(wp) :: zt, ztp, zsk      ! dummy variables
107      REAL(wp) :: zsol              ! solubility
108      REAL(wp) :: zsch              ! schmidt number
109      REAL(wp) :: zv2               ! wind speed ( square)
110      REAL(wp) :: zpv               ! piston velocity
111      REAL(wp) :: zdemi, ztra
112      REAL(wp), POINTER, DIMENSION(:,:  ) :: zatmbc14 
113      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdecay
114      !!---------------------------------------------------------------------
115      !
116      IF( nn_timing == 1 )  CALL timing_start('trc_sms_c14b')
117      !
118      ! Allocate temporary workspace
119      CALL wrk_alloc( jpi, jpj,      zatmbc14 )
120      CALL wrk_alloc( jpi, jpj, jpk, zdecay   )
121
122      IF( kt == nittrc000 )  THEN         ! Computation of decay coeffcient
123         zdemi   = 5730._wp
124         xlambda = LOG(2.) / zdemi / ( nyear_len(1) * rday )
125         xdecay  = EXP( - xlambda * rdt )
126         xaccum  = 1._wp -  xdecay
127         !
128         IF( ln_rsttr ) THEN
129            IF(lwp) WRITE(numout,*)
130            IF(lwp) WRITE(numout,*) ' Read specific variables from C14b model '
131            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
132            CALL iom_get( numrtr, jpdom_autoglo, 'qint_c14', qint_c14 )
133         ENDIF
134         !
135         IF(lwp) WRITE(numout,*)
136         !
137      ENDIF
138
139      ! Temporal interpolation
140      ! ----------------------
141      iyear_beg = nyear - 1954 + ( nyear_res_b - 1700 - nyear_beg_b  ) ! JMM Very dangerous
142                                                                       ! nyear=0001 for 1955
143      ! For Atmospheric C14 concentrations:
144      iyear_beg1 = iyear_beg
145      imonth1    = nmonth
146      IF ( imonth1 <= 6 ) THEN
147         iyear_beg1 = iyear_beg1 - 1
148         im1       =  6 - imonth1 + 1
149         im2       =  6 + imonth1 - 1
150      ELSE
151         im1       = 12 - imonth1 + 7
152         im2       =      imonth1 - 7
153      ENDIF
154      iyear_end1 = iyear_beg1 + 1
155
156      iyear_beg1 = MAX( iyear_beg1, 1   )
157      iyear_end1 = MIN( iyear_end1, 241 )
158
159
160      ! For Atmospheric CO2 concentrations:
161      iyear_beg2 = 2 * iyear_beg - 1
162      imonth2    = INT( nmonth / 2. )
163
164      IF ( imonth2 <= 3 ) THEN
165         iyear_beg2 = iyear_beg2 - 1 
166         in1       = 3 - imonth2 + 1
167         in2       = 3 + imonth2 - 1
168      ELSE
169         in1       = 6 - imonth2 + 4
170         in2       =     imonth2 - 4
171      ENDIF
172     
173      iyear_end2 = iyear_beg2 + 1
174     
175      iyear_beg2 = MAX( iyear_beg2, 1   )
176      iyear_end2 = MIN( iyear_end2, 482 )
177
178
179      !  ----------------------------------------------------------------
180      !  As explained by the TDB 89 papers, C-14/C-12 is the same
181      !  as C-14 concentration in this special case (no fractionation
182      !  effects in this model, which thus allow direct comparison
183      !  to del-C-14, after unit change from above).
184      ! -----------------------------------------------------------------------
185      !  Calc C-14 in atmosphere based on interp of IPCC (M. Heimann) data
186      !        - Compare input to calculated C-14 for each time in record
187      !-------------------------------------------------------------------------
188     
189     
190      !  Temporal and spatial interpolation at time k
191      ! --------------------------------------------------
192     
193      ! Compute atmospheric C-14 for each zone (90-20S, 20S-20N, 20-90N)
194      DO jz = 1, jpzon
195         zonbc14(jz) = (  bomb(iyear_beg1,jz) * FLOAT( im1 )  &
196              &         + bomb(iyear_end1,jz) * FLOAT( im2 ) ) / 12.
197         ! C-14 exactly starts at zero :
198         ! JMM +Zouhair : Slightly negative values are set to 0 (when perturbation approaches)
199         zonbc14(jz) = MAX( zonbc14(jz), 0. )
200      END DO
201     
202     
203      !  For each (i,j)-box, with information from the fractional area
204      !  (zonmean), computes area-weighted mean to give the atmospheric C-14
205      !  ----------------------------------------------------------------
206      zatmbc14(:,:) = zonbc14(1) * fareaz(:,:,1)   &
207         &          + zonbc14(2) * fareaz(:,:,2)   &
208         &          + zonbc14(3) * fareaz(:,:,3)
209     
210      ! time interpolation of CO2 concentrations to it time step 
211      zpco2at = (  spco2(iyear_beg2) * FLOAT( in1 )              &
212           &     + spco2(iyear_end2) * FLOAT( in2 ) ) / 6.
213
214      IF(lwp) THEN
215          WRITE(numout, *) 'time : ', kt, ' CO2 year begin/end :',iyear_beg2,'/',iyear_end2,   &
216          &                ' CO2 concen : ',zpco2at 
217          WRITE(numout, *) 'time : ', kt, ' C14 year begin/end :',iyear_beg1,'/',iyear_end1,   &
218          &                ' C14B concen (Z1/Z2/Z3) : ',zonbc14(1),'/',zonbc14(2),'/',zonbc14(3)
219      ENDIF
220
221      !  Determine seasonally variable gas exchange coefficient
222      !----------------------------------------------------------
223      !  Computes the Schmidt number of CO2 in seawater using the formulation
224      !  presented by Wanninkhof (1992, J. Geophys. Res., 97,7373-7382). 
225      ! -------------------------------------------------------------------
226      DO jj = 1, jpj
227         DO ji = 1, jpi 
228            !  Computation of solubility 
229            IF (tmask(ji,jj,1) >  0.) THEN
230               ztp  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
231               zsk  = 0.023517 + ztp * ( -0.023656 + 0.0047036 * ztp )
232               zsol = EXP( -60.2409 + 93.4517 / ztp  + 23.3585 * LOG( ztp ) + zsk * tsn(ji,jj,1,jp_sal) )
233               ! convert solubilities [mol/(l * atm)] -> [mol/(m^3 * ppm)]
234               zsol = zsol * 1.e-03
235            ELSE
236               zsol = 0._wp
237            ENDIF
238
239            ! speed transfert : Formulation of Wanninkhof (1992, JGR, 97,7373-7382)
240            ! JMM/Zouhair : coef of 0.25 rather than 0.3332 for CORe wind speed
241
242            ! Computes the Schmidt number of CO2 in seawater
243            zt   = tsn(ji,jj,1,jp_tem)
244            zsch = 2073.1 + zt * ( -125.62 + zt * (3.6276 - 0.043219 * zt ) )
245
246            ! Wanninkhof Piston velocity and convert from units [cm/hr] -> [m/s]
247            zv2  = wndm(ji,jj) * wndm(ji,jj)
248            zsch = zsch / 660.
249            zpv  = ( 0.25  * zv2 / SQRT(zsch) ) * xconv2 * tmask(ji,jj,1)
250
251            ! Flux of Bomb C-14 from air-sea : speed*(conc. at equil-conc. at surf)
252            ! in C-14 (orr-model-units) / m**2 * s
253            qtr_c14(ji,jj) = -zpv * zsol * zpco2at                                   &
254                 &                       * ( trb(ji,jj,1,jpc14) - zatmbc14(ji,jj) )  &
255#if defined key_degrad
256                 &                       * facvol(ji,jj,1)                           &
257#endif
258                  &                      * tmask(ji,jj,1) * ( 1. - fr_i(ji,jj) ) / 2.
259            ! Add the surface flux to the trend
260            tra(ji,jj,1,jpc14) = tra(ji,jj,1,jpc14) + qtr_c14(ji,jj) / fse3t(ji,jj,1) 
261           
262            ! cumulation of surface flux at each time step
263            qint_c14(ji,jj) = qint_c14(ji,jj) + qtr_c14(ji,jj) * rdt
264            !
265         END DO
266      END DO
267
268      ! Computation of decay effects on tracer concentration
269      DO jk = 1, jpk
270         DO jj = 1, jpj
271            DO ji = 1, jpi
272#if defined key_degrad
273               zdecay(ji,jj,jk) = trn(ji,jj,jk,jpc14) * ( 1. - EXP( -xlambda * rdt * facvol(ji,jj,jk) ) )
274#else
275               zdecay(ji,jj,jk) = trn(ji,jj,jk,jpc14) * xaccum
276#endif
277               tra(ji,jj,jk,jpc14) = tra(ji,jj,jk,jpc14) - zdecay(ji,jj,jk) / rdt
278               !
279            END DO
280         END DO
281      END DO
282
283      !
284      IF( lrst_trc ) THEN
285         IF(lwp) WRITE(numout,*)
286         IF(lwp) WRITE(numout,*) 'trc_sms_c14b : cumulated input function fields written in ocean restart file ',   &
287            &                    'at it= ', kt,' date= ', ndastp
288         IF(lwp) WRITE(numout,*) '~~~~'
289         CALL iom_rstput( kt, nitrst, numrtw, 'qint_c14', qint_c14 )
290      ENDIF
291      !   
292      IF( lk_iomput ) THEN
293        CALL iom_put( "qtrC14b"  , qtr_c14  )
294        CALL iom_put( "qintC14b" , qint_c14 )
295        CALL iom_put( "fdecay"   , zdecay   )
296      ELSE
297         IF( ln_diatrc ) THEN
298            trc2d(:,:  ,jp_c14b0_2d     ) = qtr_c14 (:,:)
299            trc2d(:,:  ,jp_c14b0_2d + 1 ) = qint_c14(:,:)
300            trc3d(:,:,:,jp_c14b0_3d     ) = zdecay  (:,:,:)
301         ENDIF
302      ENDIF
303
304      IF( l_trdtrc )  CALL trd_trc( tra(:,:,:,jpc14), jpc14, jptra_sms, kt )   ! save trends
305
306      CALL wrk_dealloc( jpi, jpj,      zatmbc14 )
307      CALL wrk_dealloc( jpi, jpj, jpk, zdecay   )
308      !
309      IF( nn_timing == 1 )  CALL timing_stop('trc_sms_c14b')
310      !
311   END SUBROUTINE trc_sms_c14b
312
313
314   INTEGER FUNCTION trc_sms_c14b_alloc()
315      !!----------------------------------------------------------------------
316      !!                  ***  ROUTINE trc_sms_c14b_alloc  ***
317      !!----------------------------------------------------------------------
318      ALLOCATE( fareaz  (jpi,jpj ,jpzon) ,     &
319         &      qtr_c14 (jpi,jpj)        ,     &
320         &      qint_c14(jpi,jpj)        , STAT=trc_sms_c14b_alloc )
321         !
322      IF( trc_sms_c14b_alloc /= 0 )   CALL ctl_warn('trc_sms_c14b_alloc: failed to allocate arrays')
323      !
324   END FUNCTION trc_sms_c14b_alloc
325
326#else
327   !!----------------------------------------------------------------------
328   !!   Default option                                         Dummy module
329   !!----------------------------------------------------------------------
330CONTAINS
331   SUBROUTINE trc_sms_c14b( kt )       ! Empty routine
332      WRITE(*,*) 'trc_freons: You should not have seen this print! error?', kt
333   END SUBROUTINE trc_sms_c14b
334#endif
335
336  !!======================================================================
337END MODULE trcsms_c14b
Note: See TracBrowser for help on using the repository browser.