source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90 @ 9455

Last change on this file since 9455 was 9455, checked in by cetlod, 3 years ago

Bugfix to avoid out-of-bounds in CFC model, see ticket #2073

  • Property svn:keywords set to Id
File size: 14.6 KB
Line 
1MODULE trcsms_cfc
2   !!======================================================================
3   !!                      ***  MODULE trcsms_cfc  ***
4   !! TOP : CFC main model
5   !!======================================================================
6   !! History :  OPA  !  1999-10  (JC. Dutay)  original code
7   !!  NEMO      1.0  !  2004-03  (C. Ethe) free form + modularity
8   !!            2.0  !  2007-12  (C. Ethe, G. Madec)  reorganisation
9   !!----------------------------------------------------------------------
10#if defined key_cfc
11   !!----------------------------------------------------------------------
12   !!   'key_cfc'                                               CFC tracers
13   !!----------------------------------------------------------------------
14   !!   trc_sms_cfc  :  compute and add CFC suface forcing to CFC trends
15   !!   cfc_init     :  sets constants for CFC surface forcing computation
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_cfc         ! called in ???   
28   PUBLIC   trc_sms_cfc_alloc   ! called in trcini_cfc.F90
29
30   INTEGER , PUBLIC, PARAMETER ::   jphem  =   2   ! parameter for the 2 hemispheres
31   INTEGER , PUBLIC            ::   jpyear         ! Number of years read in input data file (in trcini_cfc)
32   INTEGER , PUBLIC            ::   ndate_beg      ! initial calendar date (aammjj) for CFC
33   INTEGER , PUBLIC            ::   nyear_res      ! restoring time constant (year)
34   INTEGER , PUBLIC            ::   nyear_beg      ! initial year (aa)
35   CHARACTER(len=200), PUBLIC  ::   clnamecfc      ! Input filename of CFCs atm. concentrations
36   
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   p_cfc    ! partial hemispheric pressure for CFC
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   xphem    ! spatial interpolation factor for patm
39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qtr_cfc  ! flux at surface
40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qint_cfc ! cumulative flux
41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   atm_cfc  ! partial hemispheric pressure for used CFC
42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   patm     ! atmospheric function
43
44   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   soa      ! coefficient for solubility of CFC [mol/l/atm]
45   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sob      !    "               "
46   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sca      ! coefficients for schmidt number in degrees Celsius
47   !                          ! coefficients for conversion
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   REAL(wp) ::   xconv4 = 1.0e-12      ! conversion from mol/m3/atm to mol/m3/pptv
52
53   !! * Substitutions
54#  include "top_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
57   !! $Id$
58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE trc_sms_cfc( kt )
63      !!----------------------------------------------------------------------
64      !!                     ***  ROUTINE trc_sms_cfc  ***
65      !!
66      !! ** Purpose :   Compute the surface boundary contition on CFC 11
67      !!             passive tracer associated with air-mer fluxes and add it
68      !!             to the general trend of tracers equations.
69      !!
70      !! ** Method  : - get the atmospheric partial pressure - given in pico -
71      !!              - computation of solubility ( in 1.e-12 mol/l then in 1.e-9 mol/m3)
72      !!              - computation of transfert speed ( given in cm/hour ----> cm/s )
73      !!              - the input function is given by :
74      !!                speed * ( concentration at equilibrium - concentration at surface )
75      !!              - the input function is in pico-mol/m3/s and the
76      !!                CFC concentration in pico-mol/m3
77      !!----------------------------------------------------------------------
78      !
79      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
80      !
81      INTEGER  ::   ji, jj, jn, jl, jm, js
82      INTEGER  ::   iyear_beg, iyear_end
83      INTEGER  ::   im1, im2, ierr
84      REAL(wp) ::   ztap, zdtap       
85      REAL(wp) ::   zt1, zt2, zt3, zt4, zv2
86      REAL(wp) ::   zsol      ! solubility
87      REAL(wp) ::   zsch      ! schmidt number
88      REAL(wp) ::   zpp_cfc   ! atmospheric partial pressure of CFC
89      REAL(wp) ::   zca_cfc   ! concentration at equilibrium
90      REAL(wp) ::   zak_cfc   ! transfert coefficients
91      REAL(wp), ALLOCATABLE, DIMENSION(:,:)  ::   zpatm     ! atmospheric function
92      !!----------------------------------------------------------------------
93      !
94      !
95      IF( nn_timing == 1 )  CALL timing_start('trc_sms_cfc')
96
97      IF(lwp) WRITE(numout,*)
98      IF(lwp) WRITE(numout,*) ' trc_sms_cfc:  CFC model'
99      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
100      !
101      ALLOCATE( zpatm(jphem,jp_cfc), STAT=ierr )
102      IF( ierr > 0 ) THEN
103         CALL ctl_stop( 'trc_sms_cfc: unable to allocate zpatm array' )   ;   RETURN
104      ENDIF
105
106      IF( kt == nittrc000 )   CALL cfc_init
107
108      ! Temporal interpolation
109      ! ----------------------
110      iyear_beg = nyear - 1900
111      IF ( nmonth <= 6 ) THEN
112         iyear_beg = iyear_beg - 1
113         im1       =  6 - nmonth + 1
114         im2       =  6 + nmonth - 1
115      ELSE
116         im1       = 12 - nmonth + 7
117         im2       =      nmonth - 7
118      ENDIF
119      ! Avoid bad interpolation if starting date is =< 1900
120      if( iyear_beg .LE. 0 ) iyear_beg = 1
121      !
122      iyear_end = iyear_beg + 1
123      !                                                  !------------!
124      DO jl = 1, jp_cfc                                  !  CFC loop  !
125         !                                               !------------!
126         jn = jp_cfc0 + jl - 1
127         ! time interpolation at time kt
128         DO jm = 1, jphem
129            zpatm(jm,jl) = (  atm_cfc(iyear_beg, jm, jl) * REAL(im1, wp)  &
130               &           +  atm_cfc(iyear_end, jm, jl) * REAL(im2, wp) ) / 12.
131         END DO
132         
133         !                                                         !------------!
134         DO jj = 1, jpj                                            !  i-j loop  !
135            DO ji = 1, jpi                                         !------------!
136 
137               ! space interpolation
138               zpp_cfc  =       xphem(ji,jj)   * zpatm(1,jl)   &
139                  &     + ( 1.- xphem(ji,jj) ) * zpatm(2,jl)
140
141               ! Computation of concentration at equilibrium : in picomol/l
142               ! coefficient for solubility for CFC-11/12 in  mol/l/atm
143               IF( tmask(ji,jj,1) .GE. 0.5 ) THEN
144                  ztap  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
145                  zdtap = sob(1,jl) + ztap * ( sob(2,jl) + ztap * sob(3,jl) ) 
146                  zsol  =  EXP( soa(1,jl) + soa(2,jl) / ztap + soa(3,jl) * LOG( ztap )   &
147                     &                    + soa(4,jl) * ztap * ztap + tsn(ji,jj,1,jp_sal) * zdtap ) 
148               ELSE
149                  zsol  = 0.e0
150               ENDIF
151               ! conversion from mol/l/atm to mol/m3/atm and from mol/m3/atm to mol/m3/pptv   
152               zsol = xconv4 * xconv3 * zsol * tmask(ji,jj,1) 
153               ! concentration at equilibrium
154               zca_cfc = xconv1 * zpp_cfc * zsol * tmask(ji,jj,1)             
155 
156               ! Computation of speed transfert
157               !    Schmidt number revised in Wanninkhof (2014)
158               zt1  = tsn(ji,jj,1,jp_tem)
159               zt2  = zt1 * zt1 
160               zt3  = zt1 * zt2
161               zt4  = zt2 * zt2
162               zsch = sca(1,jl) + sca(2,jl) * zt1 + sca(3,jl) * zt2 + sca(4,jl) * zt3 + sca(5,jl) * zt4
163
164               !    speed transfert : formulae revised in Wanninkhof (2014)
165               zv2     = wndm(ji,jj) * wndm(ji,jj)
166               zsch    = zsch / 660.
167               zak_cfc = ( 0.251 * xconv2 * zv2 / SQRT(zsch) ) * tmask(ji,jj,1)
168
169               ! Input function  : speed *( conc. at equil - concen at surface )
170               ! trn in pico-mol/l idem qtr; ak in en m/a
171               qtr_cfc(ji,jj,jl) = -zak_cfc * ( trb(ji,jj,1,jn) - zca_cfc )   &
172#if defined key_degrad
173                  &                         * facvol(ji,jj,1)                           &
174#endif
175                  &                         * tmask(ji,jj,1) * ( 1. - fr_i(ji,jj) )
176               ! Add the surface flux to the trend
177               tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + qtr_cfc(ji,jj,jl) / fse3t(ji,jj,1) 
178
179               ! cumulation of surface flux at each time step
180               qint_cfc(ji,jj,jl) = qint_cfc(ji,jj,jl) + qtr_cfc(ji,jj,jl) * rdt
181               !                                               !----------------!
182            END DO                                             !  end i-j loop  !
183         END DO                                                !----------------!
184         !                                                  !----------------!
185      END DO                                                !  end CFC loop  !
186      !
187      IF( lrst_trc ) THEN
188         IF(lwp) WRITE(numout,*)
189         IF(lwp) WRITE(numout,*) 'trc_sms_cfc : cumulated input function fields written in ocean restart file ',   &
190            &                    'at it= ', kt,' date= ', ndastp
191         IF(lwp) WRITE(numout,*) '~~~~'
192         jl = 0
193         DO jn = jp_cfc0, jp_cfc1
194            jl = jl + 1
195            CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
196         END DO
197      ENDIF                                           
198      !
199      IF( lk_iomput ) THEN
200         jl = 0
201         DO jn = jp_cfc0, jp_cfc1
202            jl = jl + 1
203            CALL iom_put( 'qtr_'//TRIM(ctrcnm(jn)) , qtr_cfc (:,:,jl) )
204            CALL iom_put( 'qint_'//TRIM(ctrcnm(jn)), qint_cfc(:,:,jl) )
205         ENDDO
206      ELSE
207         IF( ln_diatrc ) THEN
208            jl = 0
209            DO jn = jp_cfc0_2d, jp_cfc1_2d, 2
210               jl = jl + 1
211               trc2d(:,:,jn    ) = qtr_cfc (:,:,jl)
212               trc2d(:,:,jn + 1) = qint_cfc(:,:,jl)
213            ENDDO
214         END IF
215      END IF
216      !
217      IF( l_trdtrc ) THEN
218          DO jn = jp_cfc0, jp_cfc1
219            CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
220          END DO
221      END IF
222      !
223      IF( nn_timing == 1 )  CALL timing_stop('trc_sms_cfc')
224      !
225   END SUBROUTINE trc_sms_cfc
226
227
228   SUBROUTINE cfc_init
229      !!---------------------------------------------------------------------
230      !!                     ***  cfc_init  *** 
231      !!
232      !! ** Purpose : sets constants for CFC model
233      !!---------------------------------------------------------------------
234      INTEGER :: jn, jl
235      !!----------------------------------------------------------------------
236      !
237      jn = 0 
238      ! coefficient for CFC11
239      !----------------------
240      if ( lp_cfc11 ) then
241         jn = jn + 1
242         ! Solubility
243         soa(1,jn) = -229.9261 
244         soa(2,jn) =  319.6552
245         soa(3,jn) =  119.4471
246         soa(4,jn) =  -1.39165
247
248         sob(1,jn) =  -0.142382
249         sob(2,jn) =   0.091459
250         sob(3,jn) =  -0.0157274
251
252         ! Schmidt number
253         sca(1,jn) = 3579.2
254         sca(2,jn) = -222.63
255         sca(3,jn) = 7.5749
256         sca(4,jn) = -0.14595
257         sca(5,jn) = 0.0011874
258
259         ! atm. concentration
260         atm_cfc(:,:,jn) = p_cfc(:,:,1)
261      endif
262
263      ! coefficient for CFC12
264      !----------------------
265      if ( lp_cfc12 ) then
266         jn = jn + 1
267         ! Solubility
268         soa(1,jn) = -218.0971
269         soa(2,jn) =  298.9702
270         soa(3,jn) =  113.8049
271         soa(4,jn) =  -1.39165
272
273         sob(1,jn) =  -0.143566
274         sob(2,jn) =   0.091015
275         sob(3,jn) =  -0.0153924
276
277         ! schmidt number
278         sca(1,jn) = 3828.1
279         sca(2,jn) = -249.86
280         sca(3,jn) = 8.7603
281         sca(4,jn) = -0.1716
282         sca(5,jn) = 0.001408
283
284         ! atm. concentration
285         atm_cfc(:,:,jn) = p_cfc(:,:,2)
286      endif
287
288      ! coefficient for SF6
289      !----------------------
290      if ( lp_sf6 ) then
291         jn = jn + 1
292         ! Solubility
293         soa(1,jn) = -80.0343
294         soa(2,jn) = 117.232
295         soa(3,jn) =  29.5817
296         soa(4,jn) =   0.0
297
298         sob(1,jn) =  0.0335183 
299         sob(2,jn) = -0.0373942 
300         sob(3,jn) =  0.00774862
301
302         ! schmidt number
303         sca(1,jn) = 3177.5
304         sca(2,jn) = -200.57
305         sca(3,jn) = 6.8865
306         sca(4,jn) = -0.13335
307         sca(5,jn) = 0.0010877
308 
309         ! atm. concentration
310         atm_cfc(:,:,jn) = p_cfc(:,:,3)
311       endif
312
313      IF( ln_rsttr ) THEN
314         IF(lwp) WRITE(numout,*)
315         IF(lwp) WRITE(numout,*) ' Read specific variables from CFC model '
316         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
317         !
318         jl = 0
319         DO jn = jp_cfc0, jp_cfc1
320            jl = jl + 1
321            CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) ) 
322         END DO
323      ENDIF
324      IF(lwp) WRITE(numout,*)
325      !
326   END SUBROUTINE cfc_init
327
328
329   INTEGER FUNCTION trc_sms_cfc_alloc()
330      !!----------------------------------------------------------------------
331      !!                     ***  ROUTINE trc_sms_cfc_alloc  ***
332      !!----------------------------------------------------------------------
333      ALLOCATE( xphem   (jpi,jpj)        , atm_cfc(jpyear,jphem,jp_cfc)  ,    &
334         &      qtr_cfc (jpi,jpj,jp_cfc) , qint_cfc(jpi,jpj,jp_cfc)      ,    &
335         &      soa(4,jp_cfc)    ,  sob(3,jp_cfc)   ,  sca(5,jp_cfc)     ,    &
336         &      STAT=trc_sms_cfc_alloc )
337         !
338      IF( trc_sms_cfc_alloc /= 0 ) CALL ctl_warn('trc_sms_cfc_alloc : failed to allocate arrays.')
339      !
340   END FUNCTION trc_sms_cfc_alloc
341
342#else
343   !!----------------------------------------------------------------------
344   !!   Dummy module                                         No CFC tracers
345   !!----------------------------------------------------------------------
346CONTAINS
347   SUBROUTINE trc_sms_cfc( kt )       ! Empty routine
348      WRITE(*,*) 'trc_sms_cfc: You should not have seen this print! error?', kt
349   END SUBROUTINE trc_sms_cfc
350#endif
351
352   !!======================================================================
353END MODULE trcsms_cfc
Note: See TracBrowser for help on using the repository browser.