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/UKMO/r5936_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/TOP_SRC/C14b – NEMO

source: branches/UKMO/r5936_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90 @ 7113

Last change on this file since 7113 was 7113, checked in by jcastill, 7 years ago

Remove again the svn keywords, as it did not work before

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