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/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/C14b – NEMO

source: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90 @ 2830

Last change on this file since 2830 was 2830, checked in by kpedwards, 13 years ago

Updates to average physics variables for TOP substepping.

  • Property svn:keywords set to Id
File size: 14.3 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 trdmod_oce
21   USE trdmod_trc
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   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/SMS/trcc14bomb.F90,v 1.2 2005/11/14 16:42:43 opalod Exp $
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      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
97      USE wrk_nemo, ONLY:   zatmbc14 => wrk_2d_1
98      USE wrk_nemo, ONLY:   zw3d     => wrk_3d_1
99      !
100      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
101      !
102      INTEGER :: ji, jj, jk, jz     ! dummy loop indices
103      INTEGER :: iyear_beg , iyear_beg1, iyear_end1 
104      INTEGER :: iyear_beg2, iyear_end2 
105      INTEGER :: imonth1, im1, in1 
106      INTEGER :: imonth2, im2, in2 
107      REAL(wp), DIMENSION(jpzon) :: zonbc14       !: time interp atm C14
108      REAL(wp)                   :: zpco2at       !: time interp atm C02
109      REAL(wp) :: zt, ztp, zsk      ! dummy variables
110      REAL(wp) :: zsol              ! solubility
111      REAL(wp) :: zsch              ! schmidt number
112      REAL(wp) :: zv2               ! wind speed ( square)
113      REAL(wp) :: zpv               ! piston velocity
114      REAL(wp) :: zdemi, ztra
115      !!----------------------------------------------------------------------
116
117      IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 1) ) THEN
118         CALL ctl_stop('trc_sms_c14b : requested workspace arrays unavailable')   ;   RETURN
119      ENDIF
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      ENDIF
127
128      ! Temporal interpolation
129      ! ----------------------
130      iyear_beg = nyear - 1954 + ( nyear_res_b - 1700 - nyear_beg_b  ) ! JMM Very dangerous
131                                                                       ! nyear=0001 for 1955
132      ! For Atmospheric C14 concentrations:
133      iyear_beg1 = iyear_beg
134      imonth1    = nmonth
135      IF ( imonth1 <= 6 ) THEN
136         iyear_beg1 = iyear_beg1 - 1
137         im1       =  6 - imonth1 + 1
138         im2       =  6 + imonth1 - 1
139      ELSE
140         im1       = 12 - imonth1 + 7
141         im2       =      imonth1 - 7
142      ENDIF
143      iyear_end1 = iyear_beg1 + 1
144
145      iyear_beg1 = MAX( iyear_beg1, 1   )
146      iyear_end1 = MIN( iyear_end1, 241 )
147
148
149      ! For Atmospheric CO2 concentrations:
150      iyear_beg2 = 2 * iyear_beg - 1
151      imonth2    = INT( nmonth / 2. )
152
153      IF ( imonth2 <= 3 ) THEN
154         iyear_beg2 = iyear_beg2 - 1 
155         in1       = 3 - imonth2 + 1
156         in2       = 3 + imonth2 - 1
157      ELSE
158         in1       = 6 - imonth2 + 4
159         in2       =     imonth2 - 4
160      ENDIF
161     
162      iyear_end2 = iyear_beg2 + 1
163     
164      iyear_beg2 = MAX( iyear_beg2, 1   )
165      iyear_end2 = MIN( iyear_end2, 482 )
166
167
168      !  ----------------------------------------------------------------
169      !  As explained by the TDB 89 papers, C-14/C-12 is the same
170      !  as C-14 concentration in this special case (no fractionation
171      !  effects in this model, which thus allow direct comparison
172      !  to del-C-14, after unit change from above).
173      ! -----------------------------------------------------------------------
174      !  Calc C-14 in atmosphere based on interp of IPCC (M. Heimann) data
175      !        - Compare input to calculated C-14 for each time in record
176      !-------------------------------------------------------------------------
177     
178     
179      !  Temporal and spatial interpolation at time k
180      ! --------------------------------------------------
181     
182      ! Compute atmospheric C-14 for each zone (90-20S, 20S-20N, 20-90N)
183      DO jz = 1, jpzon
184         zonbc14(jz) = (  bomb(iyear_beg1,jz) * FLOAT( im1 )  &
185              &         + bomb(iyear_end1,jz) * FLOAT( im2 ) ) / 12.
186         ! C-14 exactly starts at zero :
187         ! JMM +Zouhair : Slightly negative values are set to 0 (when perturbation approaches)
188         zonbc14(jz) = MAX( zonbc14(jz), 0. )
189      END DO
190     
191     
192      !  For each (i,j)-box, with information from the fractional area
193      !  (zonmean), computes area-weighted mean to give the atmospheric C-14
194      !  ----------------------------------------------------------------
195      zatmbc14(:,:) = zonbc14(1) * fareaz(:,:,1)   &
196         &          + zonbc14(2) * fareaz(:,:,2)   &
197         &          + zonbc14(3) * fareaz(:,:,3)
198     
199      ! time interpolation of CO2 concentrations to it time step 
200      zpco2at = (  spco2(iyear_beg2) * FLOAT( in1 )              &
201           &     + spco2(iyear_end2) * FLOAT( in2 ) ) / 6.
202
203      IF(lwp) THEN
204          WRITE(numout, *) 'time : ', kt, ' CO2 year begin/end :',iyear_beg2,'/',iyear_end2,   &
205          &                ' CO2 concen : ',zpco2at 
206          WRITE(numout, *) 'time : ', kt, ' C14 year begin/end :',iyear_beg1,'/',iyear_end1,   &
207          &                ' C14B concen (Z1/Z2/Z3) : ',zonbc14(1),'/',zonbc14(2),'/',zonbc14(3)
208      ENDIF
209
210      !  Determine seasonally variable gas exchange coefficient
211      !----------------------------------------------------------
212      !  Computes the Schmidt number of CO2 in seawater using the formulation
213      !  presented by Wanninkhof (1992, J. Geophys. Res., 97,7373-7382). 
214      ! -------------------------------------------------------------------
215      DO jj = 1, jpj
216         DO ji = 1, jpi 
217            !  Computation of solubility 
218            IF (tmask(ji,jj,1) >  0.) THEN
219               ztp  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
220               zsk  = 0.023517 + ztp * ( -0.023656 + 0.0047036 * ztp )
221               zsol = EXP( -60.2409 + 93.4517 / ztp  + 23.3585 * LOG( ztp ) + zsk * tsn(ji,jj,1,jp_sal) )
222               ! convert solubilities [mol/(l * atm)] -> [mol/(m^3 * ppm)]
223               zsol = zsol * 1.e-03
224            ELSE
225               zsol = 0._wp
226            ENDIF
227
228            ! speed transfert : Formulation of Wanninkhof (1992, JGR, 97,7373-7382)
229            ! JMM/Zouhair : coef of 0.25 rather than 0.3332 for CORe wind speed
230
231            ! Computes the Schmidt number of CO2 in seawater
232            zt   = tsn(ji,jj,1,jp_tem)
233            zsch = 2073.1 + zt * ( -125.62 + zt * (3.6276 - 0.043219 * zt ) )
234
235            ! Wanninkhof Piston velocity and convert from units [cm/hr] -> [m/s]
236            zv2  = wndm(ji,jj) * wndm(ji,jj)
237            zsch = zsch / 660.
238            zpv  = ( 0.25  * zv2 / SQRT(zsch) ) * xconv2 * tmask(ji,jj,1)
239
240            ! Flux of Bomb C-14 from air-sea : speed*(conc. at equil-conc. at surf)
241            ! in C-14 (orr-model-units) / m**2 * s
242            qtr_c14(ji,jj) = -zpv * zsol * zpco2at                                   &
243                 &                       * ( trb(ji,jj,1,jpc14) - zatmbc14(ji,jj) )  &
244#if defined key_degrad
245                 &                       * facvol(ji,jj,1)                           &
246#endif
247                  &                      * tmask(ji,jj,1) * ( 1. - fr_i(ji,jj) ) / 2.
248
249            ! Add the surface flux to the trend
250            tra(ji,jj,1,jpc14) = tra(ji,jj,1,jpc14) + qtr_c14(ji,jj) / fse3t(ji,jj,1) 
251           
252            ! cumulation of surface flux at each time step
253            qint_c14(ji,jj) = qint_c14(ji,jj) + qtr_c14(ji,jj) * rdt
254
255# if defined key_diatrc && ! defined key_iomput
256            ! Save 2D diagnostics
257            trc2d(ji,jj,jp_c14b0_2d    ) = qtr_c14 (ji,jj)
258            trc2d(ji,jj,jp_c14b0_2d + 1) = qint_c14(ji,jj)
259# endif
260         END DO
261      END DO
262
263      ! Computation of decay effects on tracer concentration
264      DO jk = 1, jpk
265         DO jj = 1, jpj
266            DO ji = 1, jpi
267#if ! defined key_degrad
268               ztra = trn(ji,jj,jk,jpc14) * xaccum
269#else
270               ztra = trn(ji,jj,jk,jpc14) * ( 1. - EXP( -xlambda * rdt * facvol(ji,jj,jk) ) )
271#endif
272               tra(ji,jj,jk,jpc14) = tra(ji,jj,jk,jpc14) - ztra / rdt
273#if defined key_diatrc
274               ! Save 3D diagnostics
275# if ! defined key_iomput
276               trc3d(ji,jj,jk,jp_c14b0_3d ) = ztra    !  radioactive decay
277# else
278               zw3d(ji,jj,jk) = ztra    !  radioactive decay
279# endif
280#endif
281            END DO
282         END DO
283      END DO
284
285#if defined key_diatrc  && defined key_iomput
286      CALL iom_put( "qtrC14b"  , qtr_c14  )
287      CALL iom_put( "qintC14b" , qint_c14 )
288#endif
289#if defined key_diatrc  && defined key_iomput
290      CALL iom_put( "fdecay" , zw3d )
291#endif
292      IF( l_trdtrc )   CALL trd_mod_trc( tra(:,:,:,jpc14), jpc14, jptra_trd_sms, kt )   ! save trends
293
294      IF( wrk_not_released(2, 1) .OR.   &
295          wrk_not_released(3, 1) )   CALL ctl_stop('trc_sms_c14b : failed to release workspace arrays')
296      !
297   END SUBROUTINE trc_sms_c14b
298
299
300   INTEGER FUNCTION trc_sms_c14b_alloc()
301      !!----------------------------------------------------------------------
302      !!                  ***  ROUTINE trc_sms_c14b_alloc  ***
303      !!----------------------------------------------------------------------
304      ALLOCATE( fareaz  (jpi,jpj ,jpzon) ,     &
305         &      qtr_c14 (jpi,jpj)        ,     &
306         &      qint_c14(jpi,jpj)        , STAT=trc_sms_c14b_alloc )
307         !
308      IF( trc_sms_c14b_alloc /= 0 )   CALL ctl_warn('trc_sms_c14b_alloc: failed to allocate arrays')
309      !
310   END FUNCTION trc_sms_c14b_alloc
311
312#else
313   !!----------------------------------------------------------------------
314   !!   Default option                                         Dummy module
315   !!----------------------------------------------------------------------
316CONTAINS
317   SUBROUTINE trc_sms_c14b( kt )       ! Empty routine
318      WRITE(*,*) 'trc_freons: You should not have seen this print! error?', kt
319   END SUBROUTINE trc_sms_c14b
320#endif
321
322  !!======================================================================
323END MODULE trcsms_c14b
Note: See TracBrowser for help on using the repository browser.