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.
trcsms_c14b.F90 in branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/TOP_SRC/C14b – NEMO

source: branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90 @ 6225

Last change on this file since 6225 was 6225, checked in by jamesharle, 8 years ago

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

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