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 tags/nemo_v3_1_beta/NEMO/TOP_SRC/C14b – NEMO

source: tags/nemo_v3_1_beta/NEMO/TOP_SRC/C14b/trcsms_c14b.F90 @ 9247

Last change on this file since 9247 was 1252, checked in by cetlod, 16 years ago

add new modules for C14 bomb tracer model, see ticket:298

File size: 12.9 KB
Line 
1MODULE trcsms_c14b
2   !!======================================================================
3   !!                      ***  MODULE trcsms_c14b  ***
4   !! TOP : Bomb C14 main module
5   !!======================================================================
6   !! History     -   ! 1994-05 ( J. Orr ) origial code
7   !!            1.0  ! 2006-02 ( J.M. Molines )  Free form + modularity
8   !!            2.0  ! 2008-12 ( C. Ethe ) reorganisation
9   !!----------------------------------------------------------------------
10#if defined key_c14b
11   !!----------------------------------------------------------------------
12   !!   'key_c14b'                                         Bomb C14 tracer
13   !!----------------------------------------------------------------------
14   !!   trc_sms_c14b  :  compute and add C14 suface forcing to C14 trends
15   !!----------------------------------------------------------------------
16   USE daymod       ! calendar
17   USE oce_trc      ! Ocean variables
18   USE par_trc      ! TOP parameters
19   USE trc          ! TOP variables
20   USE trdmld_trc_oce
21   USE trdmld_trc
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC   trc_sms_c14b       ! called in ???   
28
29   !! * Module variables
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, DIMENSION(jpi,jpj ,jpzon) ::  fareaz !: Spatial Interpolation Factors
41   REAL(wp), PUBLIC, DIMENSION(jpmaxrec2)      ::  spco2  !: Atmospheric CO2
42 
43   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   qtr_c14      !: flux at surface
44   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   qint_c14     !: cumulative flux
45
46   REAL(wp) :: xlambda, xdecay, xaccum       ! C14 decay coef. 
47
48   REAL(wp) ::   xconv1 = 1.0          ! conversion from to
49   REAL(wp) ::   xconv2 = 0.01/3600.   ! conversion from cm/h to m/s:
50   REAL(wp) ::   xconv3 = 1.0e+3       ! conversion from mol/l/atm to mol/m3/atm
51
52  !! * Substitutions
53#  include "top_substitute.h90"
54
55  !!----------------------------------------------------------------------
56  !!  TOP 1.0 , LOCEAN-IPSL (2005)
57  !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/SMS/trcc14bomb.F90,v 1.2 2005/11/14 16:42:43 opalod Exp $
58  !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
59  !!----------------------------------------------------------------------
60
61CONTAINS
62
63  SUBROUTINE trc_sms_c14b( kt )
64     !!----------------------------------------------------------------------
65     !!                  ***  ROUTINE trc_sms_c14b  ***
66     !!
67     !! ** Purpose :   Compute the surface boundary contition on C14bomb
68     !!      passive tracer associated with air-mer fluxes and add it to
69     !!      the general trend of tracers equations.
70     !!
71     !! ** Original comments from J. Orr :
72     !!
73     !!      Calculates the input of Bomb C-14 to the surface layer of OPA
74     !!
75     !!      James Orr, LMCE, 28 October 1992
76     !!
77     !!      Initial approach mimics that of Toggweiler, Dixon, & Bryan (1989)
78     !!      (hereafter referred to as TDB) with constant gas exchange,
79     !!      although in this case, a perturbation approach is used for
80     !!      bomb-C14 so that both the ocean and atmosphere begin at zero.
81     !!      This saves tremendous amounts of computer time since no
82     !!      equilibrum run is first required (i.e., for natural C-14).
83     !!      Note: Many sensitivity tests can be run with this approach and
84     !!            one never has to make a run for natural C-14; otherwise,
85     !!            a run for natural C-14 must be run each time that one
86     !!            changes a model parameter!
87     !!
88     !!
89     !!      19 August 1993: Modified to consider Atmospheric C-14 fom IPCC.
90     !!      That is, the IPCC has provided a C-14 atmospheric record (courtesy
91     !!      of Martin Heimann) for model calibration.  This model spans from
92     !!      preindustrial times to present, in a format different than that
93     !!      given by TDB.  It must be converted to the ORR C-14 units used
94     !!      here, although in this case, the perturbation includes not only
95     !!      bomb C-14 but changes due to the Suess effect.
96     !!
97     !!----------------------------------------------------------------------
98     !! * Arguments
99     INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
100
101     !! * Local declarations
102     INTEGER ::  &
103       ji, jj, jk, jz
104
105     INTEGER ::   &
106       iyear_beg,              &
107       iyear_beg1, iyear_end1, &
108       imonth1, im1, im2,      &
109       iyear_beg2, iyear_end2, &
110       imonth2, in1, in2
111         
112     REAL(wp), DIMENSION(jpi,jpj) ::  & 
113       zatmbc14
114
115     REAL(wp), DIMENSION(jpzon)   ::  &
116       zonbc14              !: time interp atm C14
117
118     REAL(wp) ::     &
119       zpco2at              !: time interp atm C02
120
121     REAL(wp) ::     &      !: dummy variables
122       zt, ztp, zsk, &
123       zsol ,        &      !: solubility
124       zsch ,        &      !: schmidt number
125       zv2  ,        &      !: wind speed ( square)
126       zpv  ,        &      !: piston velocity
127       zdemi, ztra
128      !!----------------------------------------------------------------------
129
130      IF( kt == nittrc000 )  THEN
131         ! Computation of decay coeffcient
132         zdemi   = 5730.
133         xlambda = LOG(2.) / zdemi / raass
134         xdecay  = EXP( - xlambda * rdt )
135         xaccum  = 1.0 -  xdecay
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      DO jj = 1, jpj
206         DO ji = 1, jpi
207            zatmbc14(ji,jj) =   zonbc14(1) * fareaz(ji,jj,1)  &
208                 &           +  zonbc14(2) * fareaz(ji,jj,2)  &
209                 &           +  zonbc14(3) * fareaz(ji,jj,3)
210         END DO
211      END DO
212     
213      ! time interpolation of CO2 concentrations to it time step 
214      zpco2at = (  spco2(iyear_beg2) * FLOAT( in1 )              &
215           &     + spco2(iyear_end2) * FLOAT( in2 ) ) / 6.
216
217      IF (lwp) THEN
218          WRITE(numout, *) 'time : ', kt, ' CO2 year begin/end :',iyear_beg2,'/',iyear_end2,   &
219          &                ' CO2 concen : ',zpco2at 
220          WRITE(numout, *) 'time : ', kt, ' C14 year begin/end :',iyear_beg1,'/',iyear_end1,   &
221          &                ' C14B concen (Z1/Z2/Z3) : ',zonbc14(1),'/',zonbc14(2),'/',zonbc14(3)
222      ENDIF
223
224      !  Determine seasonally variable gas exchange coefficient
225      !----------------------------------------------------------
226      !  Computes the Schmidt number of CO2 in seawater using the formulation
227      !  presented by Wanninkhof (1992, J. Geophys. Res., 97,7373-7382). 
228      ! -------------------------------------------------------------------
229      DO jj = 1, jpj
230         DO ji = 1, jpi 
231            !  Computation of solubility 
232            IF (tmask(ji,jj,1) >  0.) THEN
233               ztp  = ( tn(ji,jj,1) + 273.16 ) * 0.01
234               zsk  = 0.023517 + ztp * ( -0.023656 + 0.0047036 * ztp )
235               zsol = EXP( -60.2409 + 93.4517 / ztp  + 23.3585 * LOG( ztp ) + zsk * sn(ji,jj,1) )
236               ! convert solubilities [mol/(l * atm)] -> [mol/(m^3 * ppm)]
237               zsol = zsol * 1.0e-03
238            ELSE
239               zsol = 0.
240            ENDIF
241
242            ! speed transfert : Formulation of Wanninkhof (1992, JGR, 97,7373-7382)
243            ! JMM/Zouhair : coef of 0.25 rather than 0.3332 for CORe wind speed
244
245            ! Computes the Schmidt number of CO2 in seawater
246            zt   = tn(ji,jj,1)
247            zsch = 2073.1 + zt * ( -125.62 + zt * (3.6276 - 0.043219 * zt ) )
248
249            ! Wanninkhof Piston velocity and convert from units [cm/hr] -> [m/s]
250            zv2  = wndm(ji,jj) * wndm(ji,jj)
251            zsch = zsch / 660.
252            zpv  = ( 0.25  * zv2 / SQRT(zsch) ) * xconv2 * tmask(ji,jj,1)
253
254            ! Flux of Bomb C-14 from air-sea : speed*(conc. at equil-conc. at surf)
255            ! in C-14 (orr-model-units) / m**2 * s
256            qtr_c14(ji,jj) = -zpv * zsol * zpco2at                                   &
257                 &                       * ( trb(ji,jj,1,jpc14) - zatmbc14(ji,jj) )  &
258#if defined key_off_degrad
259                 &                       * facvol(ji,jj,1)                           &
260#endif
261                  &                      * tmask(ji,jj,1) * ( 1. - fr_i(ji,jj) ) / 2.
262
263            ! Add the surface flux to the trend
264            tra(ji,jj,1,jpc14) = tra(ji,jj,1,jpc14) + qtr_c14(ji,jj) / fse3t(ji,jj,1) 
265           
266            ! cumulation of surface flux at each time step
267            qint_c14(ji,jj) = qint_c14(ji,jj) + qtr_c14(ji,jj) * rdt
268
269# if defined key_trc_diaadd
270            ! Save 2D diagnostics
271            trc2d(ji,jj,jp_c14b0_2d    ) = qtr_c14 (ji,jj)
272            trc2d(ji,jj,jp_c14b0_2d + 1) = qint_c14(ji,jj)
273# endif
274         END DO
275      END DO
276
277      ! Computation of decay effects on tracer concentration
278      DO jk = 1, jpk
279         DO jj = 1, jpj
280            DO ji = 1, jpi
281#if ! defined key_off_degrad
282               ztra = trn(ji,jj,jk,jpc14) * xaccum
283#else
284               ztra = trn(ji,jj,jk,jpc14) * ( 1. - EXP( -xlambda * rdt * facvol(ji,jj,jk) ) )
285#endif
286               tra(ji,jj,jk,jpc14) = tra(ji,jj,jk,jpc14) - ztra / rdt
287
288# if defined key_trc_dia3d
289               ! Save 3D diagnostics
290               trc3d(ji,jj,jk,jp_c14b0_3d ) = ztra    !  radioactive decay
291# endif
292            END DO
293         END DO
294      END DO
295
296      IF( l_trdtrc ) THEN
297         CALL trd_mod_trc( tra(:,:,:,jpc14), jpc14, jptrc_trd_sms, kt )   ! save trends
298      END IF
299
300    END SUBROUTINE trc_sms_c14b
301
302#else
303    !!----------------------------------------------------------------------
304    !!   Default option                                         Dummy module
305    !!----------------------------------------------------------------------
306CONTAINS
307  SUBROUTINE trc_sms_c14b( kt )       ! Empty routine
308    WRITE(*,*) 'trc_freons: You should not have seen this print! error?', kt
309  END SUBROUTINE trc_sms_c14b
310#endif
311
312  !!======================================================================
313END MODULE trcsms_c14b
Note: See TracBrowser for help on using the repository browser.