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_cfc.F90 in branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/CFC – NEMO

source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90 @ 8075

Last change on this file since 8075 was 8075, checked in by jpalmier, 8 years ago

JPALM -- update CFCs - add SF6 and update gas transfert param

File size: 19.0 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   !!                 !  2016-06  (J. Palmieri)  update for UKESM1
10   !!                 !  2017-04  (A. Yool)  update to add SF6, fix coefficients
11   !!----------------------------------------------------------------------
12#if defined key_cfc
13   !!----------------------------------------------------------------------
14   !!   'key_cfc'                                               CFC tracers
15   !!----------------------------------------------------------------------
16   !!   trc_sms_cfc  :  compute and add CFC suface forcing to CFC trends
17   !!   cfc_init     :  sets constants for CFC surface forcing computation
18   !!----------------------------------------------------------------------
19   USE dom_oce       ! ocean space and time domain
20   USE oce_trc       ! Ocean variables
21   USE par_trc       ! TOP parameters
22   USE trc           ! TOP variables
23   USE trd_oce
24   USE trdtrc
25   USE iom           ! I/O library
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   trc_sms_cfc         ! called in ???   
31   PUBLIC   trc_sms_cfc_alloc   ! called in trcini_cfc.F90
32
33   INTEGER , PUBLIC, PARAMETER ::   jphem  =   2   ! parameter for the 2 hemispheres
34   INTEGER , PUBLIC            ::   jpyear         ! Number of years read in CFC1112 file
35   INTEGER , PUBLIC            ::   ndate_beg      ! initial calendar date (aammjj) for CFC
36   INTEGER , PUBLIC            ::   simu_type      ! Kind of simulation: 1- Spin-up
37                                                   !                     2- Hindcast/projection
38   INTEGER , PUBLIC            ::   nyear_res      ! restoring time constant (year)
39   INTEGER , PUBLIC            ::   nyear_beg      ! initial year (aa)
40   
41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   p_cfc    ! partial hemispheric pressure for CFC
42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   xphem    ! spatial interpolation factor for patm
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qtr_cfc  ! flux at surface
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qint_cfc ! cumulative flux
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   patm     ! atmospheric function
46
47   REAL(wp), DIMENSION(4,3) ::   soa   ! coefficient for solubility of CFC [mol/l/atm]
48   REAL(wp), DIMENSION(3,3) ::   sob   !    "               "
49   REAL(wp), DIMENSION(5,3) ::   sca   ! coefficients for schmidt number in degre Celcius
50     
51   !                          ! coefficients for conversion
52   REAL(wp) ::   xconv1 = 1.0          ! conversion from to
53   REAL(wp) ::   xconv2 = 0.01/3600.   ! conversion from cm/h to m/s:
54   REAL(wp) ::   xconv3 = 1.0e+3       ! conversion from mol/l/atm to mol/m3/atm
55   REAL(wp) ::   xconv4 = 1.0e-12      ! conversion from mol/m3/atm to mol/m3/pptv
56
57   !! * Substitutions
58#  include "top_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE trc_sms_cfc( kt )
67      !!----------------------------------------------------------------------
68      !!                     ***  ROUTINE trc_sms_cfc  ***
69      !!
70      !! ** Purpose :   Compute the surface boundary contition on CFC 11
71      !!             passive tracer associated with air-mer fluxes and add it
72      !!             to the general trend of tracers equations.
73      !!
74      !! ** Method  : - get the atmospheric partial pressure - given in pico -
75      !!              - computation of solubility ( in 1.e-12 mol/l then in 1.e-9 mol/m3)
76      !!              - computation of transfert speed ( given in cm/hour ----> cm/s )
77      !!              - the input function is given by :
78      !!                speed * ( concentration at equilibrium - concentration at surface )
79      !!              - the input function is in pico-mol/m3/s and the
80      !!                CFC concentration in pico-mol/m3
81      !!----------------------------------------------------------------------
82      !
83      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
84      !
85      INTEGER  ::   ji, jj, jn, jl, jm, js
86      INTEGER  ::   iyear_beg, iyear_end, iyear_tmp
87      INTEGER  ::   im1, im2, ierr
88      REAL(wp) ::   ztap, zdtap       
89      REAL(wp) ::   zt1, zt2, zt3, zt4, zv2
90      REAL(wp) ::   zsol      ! solubility
91      REAL(wp) ::   zsch      ! schmidt number
92      REAL(wp) ::   zpp_cfc   ! atmospheric partial pressure of CFC
93      REAL(wp) ::   zca_cfc   ! concentration at equilibrium
94      REAL(wp) ::   zak_cfc   ! transfert coefficients
95      REAL(wp), ALLOCATABLE, DIMENSION(:,:)  ::   zpatm     ! atmospheric function
96      !!----------------------------------------------------------------------
97      !
98      !
99      IF( nn_timing == 1 )  CALL timing_start('trc_sms_cfc')
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      !! JPALM -- 15-06-2016 -- define 2 kinds of CFC run:
111      !!                     1- the SPIN-UP and 2- Hindcast/Projections
112      !!                     -- main difference is the way to define the year of
113      !!                     simulation, that determine the atm pCFC.
114      !!                     1-- Spin-up: our atm forcing is of 30y we cycle on.
115      !!                     So we do 90y CFC cycles to be in good
116      !!                     correspondence with the atmosphere
117      !!                     2-- Hindcast/proj, instead of nyear-1900 we keep
118      !!                     the 2 last digit, and enable 3 cycle from 1800 to 2100. 
119      !!----------------------------------------------------------------------
120      IF (simu_type==1) THEN
121         !! 1 -- SPIN-UP
122         iyear_tmp = nyear - nyear_res  !! JPALM -- in our spin-up, nyear_res is 1000
123         iyear_beg = MOD( iyear_tmp , 90 )
124         !! JPALM -- the pCFC file only got 78 years.
125         !!       So if iyear_beg > 78 then we set pCFC to 0
126         !!             iyear_beg = 0 as well -- must try to avoid obvious problems
127         !!             as Pcfc is set to 0.00 up to year 32, let set iyear_beg to year 10
128         !!          else, must add 30 to iyear_beg to match with P_cfc indices
129         !!---------------------------------------
130         IF ((iyear_beg > 77) .OR. (iyear_beg==0)) THEN
131            iyear_beg = 10
132         ELSE
133            iyear_beg = iyear_beg + 30
134         ENDIF
135      ELSEIF (simu_type==2) THEN
136         !! 2 -- Hindcast/proj
137         iyear_beg = MOD(nyear, 100)
138         IF (iyear_beg < 9)  iyear_beg = iyear_beg + 100
139         !! JPALM -- Same than previously, if iyear_beg is out of P_cfc range,
140         !!       we want to set p_CFC to 0.00 --> set iyear_beg = 10
141         IF ((iyear_beg < 30) .OR. (iyear_beg > 107)) iyear_beg = 10             
142      ENDIF
143      !!
144      IF ( nmonth <= 6 ) THEN
145         iyear_beg = iyear_beg - 1
146         im1       =  6 - nmonth + 1
147         im2       =  6 + nmonth - 1
148      ELSE
149         im1       = 12 - nmonth + 7
150         im2       =      nmonth - 7
151      ENDIF
152      iyear_end = iyear_beg + 1
153
154      !                                                  !------------!
155      DO jl = 1, jp_cfc                                  !  CFC loop  !
156         !                                               !------------!
157         jn = jp_cfc0 + jl - 1
158         ! time interpolation at time kt
159         DO jm = 1, jphem
160            zpatm(jm,jl) = (  p_cfc(iyear_beg, jm, jl) * FLOAT (im1)  &
161               &           +  p_cfc(iyear_end, jm, jl) * FLOAT (im2) ) / 12.
162         END DO
163         
164         !                                                         !------------!
165         DO jj = 1, jpj                                            !  i-j loop  !
166            DO ji = 1, jpi                                         !------------!
167 
168               ! space interpolation
169               zpp_cfc  =       xphem(ji,jj)   * zpatm(1,jl)   &
170                  &     + ( 1.- xphem(ji,jj) ) * zpatm(2,jl)
171
172               ! Computation of concentration at equilibrium : in picomol/l
173               ! coefficient for solubility for CFC-11/12 in  mol/l/atm
174               IF( tmask(ji,jj,1) .GE. 0.5 ) THEN
175                  ztap  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
176                  zdtap = sob(1,jl) + ztap * ( sob(2,jl) + ztap * sob(3,jl) ) 
177                  zsol  =  EXP( soa(1,jl) + soa(2,jl) / ztap + soa(3,jl) * LOG( ztap )   &
178                     &                    + soa(4,jl) * ztap * ztap + tsn(ji,jj,1,jp_sal) * zdtap ) 
179               ELSE
180                  zsol  = 0.e0
181               ENDIF
182               ! conversion from mol/l/atm to mol/m3/atm and from mol/m3/atm to mol/m3/pptv   
183               zsol = xconv4 * xconv3 * zsol * tmask(ji,jj,1) 
184               ! concentration at equilibrium
185               zca_cfc = xconv1 * zpp_cfc * zsol * tmask(ji,jj,1)             
186 
187               ! Computation of speed transfert
188               !    Schmidt number
189               zt1  = tsn(ji,jj,1,jp_tem)
190               zt2  = zt1 * zt1 
191               zt3  = zt1 * zt2
192               zt4  = zt1 * zt3
193               zsch = sca(1,jl) + sca(2,jl) * zt1 + sca(3,jl) * zt2 + sca(4,jl) * zt3 + sca(5,jl) * zt4
194
195               !    speed transfert : formulae of wanninkhof 1992
196               zv2     = wndm(ji,jj) * wndm(ji,jj)
197               zsch    = zsch / 660.
198               ! AXY (25/04/17): OMIP protocol specifies lower Wanninkhof (2014) value
199               ! zak_cfc = ( 0.39 * xconv2 * zv2 / SQRT(zsch) ) * tmask(ji,jj,1)
200               zak_cfc = ( 0.251 * xconv2 * zv2 / SQRT(zsch) ) * tmask(ji,jj,1)
201
202               ! Input function  : speed *( conc. at equil - concen at surface )
203               ! trn in pico-mol/l idem qtr; ak in en m/a
204               qtr_cfc(ji,jj,jl) = -zak_cfc * ( trb(ji,jj,1,jn) - zca_cfc )   &
205#if defined key_degrad
206                  &                         * facvol(ji,jj,1)                           &
207#endif
208                  &                         * tmask(ji,jj,1) * ( 1. - fr_i(ji,jj) )
209               ! Add the surface flux to the trend
210               tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + qtr_cfc(ji,jj,jl) / fse3t(ji,jj,1) 
211
212               ! cumulation of surface flux at each time step
213               qint_cfc(ji,jj,jl) = qint_cfc(ji,jj,jl) + qtr_cfc(ji,jj,jl) * rdt
214               !                                               !----------------!
215            END DO                                             !  end i-j loop  !
216         END DO                                                !----------------!
217         !                                                  !----------------!
218      END DO                                                !  end CFC loop  !
219         !
220      IF( kt == nittrc000 ) THEN
221         DO jl = 1, jp_cfc   
222             WRITE(NUMOUT,*) ' '
223             WRITE(NUMOUT,*) 'CFC interpolation verification '  !! Jpalm 
224             WRITE(NUMOUT,*) '################################## '
225             WRITE(NUMOUT,*) ' '
226               if (jl.EQ.1) then
227                   WRITE(NUMOUT,*) 'Traceur = CFC11: '
228               elseif (jl.EQ.2) then
229                   WRITE(NUMOUT,*) 'Traceur = CFC12: '
230               elseif (jl.EQ.3) then
231                   WRITE(NUMOUT,*) 'Traceur = SF6: '
232               endif
233             WRITE(NUMOUT,*) 'nyear    = ', nyear
234             WRITE(NUMOUT,*) 'nmonth   = ', nmonth
235             WRITE(NUMOUT,*) 'iyear_beg= ', iyear_beg
236             WRITE(NUMOUT,*) 'iyear_end= ', iyear_end
237             WRITE(NUMOUT,*) 'p_cfc(iyear_beg)= ',p_cfc(iyear_beg, 1, jl)
238             WRITE(NUMOUT,*) 'p_cfc(iyear_end)= ',p_cfc(iyear_end, 1, jl)
239             WRITE(NUMOUT,*) 'Im1= ',im1
240             WRITE(NUMOUT,*) 'Im2= ',im2
241             WRITE(NUMOUT,*) 'zpp_cfc = ',zpp_cfc
242             WRITE(NUMOUT,*) ' '
243         END DO 
244# if defined key_debug_medusa
245         CALL flush(numout)
246# endif
247      ENDIF
248        !
249      !IF( lrst_trc ) THEN
250      !   IF(lwp) WRITE(numout,*)
251      !   IF(lwp) WRITE(numout,*) 'trc_sms_cfc : cumulated input function fields written in ocean restart file ',   &
252      !      &                    'at it= ', kt,' date= ', ndastp
253      !   IF(lwp) WRITE(numout,*) '~~~~'
254      !   DO jn = jp_cfc0, jp_cfc1
255      !      CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
256      !   END DO
257      !ENDIF                                           
258      !
259      IF( lk_iomput ) THEN
260         CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
261         CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
262         CALL iom_put( "qtrCFC12"  , qtr_cfc (:,:,2) )
263         CALL iom_put( "qintCFC12" , qint_cfc(:,:,2) )
264         CALL iom_put( "qtrSF6"    , qtr_cfc (:,:,3) )
265         CALL iom_put( "qintSF6"   , qint_cfc(:,:,4) )
266      ELSE
267         IF( ln_diatrc ) THEN
268            trc2d(:,:,jp_cfc0_2d    ) = qtr_cfc (:,:,1)
269            trc2d(:,:,jp_cfc0_2d + 1) = qint_cfc(:,:,1)
270            trc2d(:,:,jp_cfc0_2d + 2) = qtr_cfc (:,:,2)
271            trc2d(:,:,jp_cfc0_2d + 3) = qint_cfc(:,:,2)
272            trc2d(:,:,jp_cfc0_2d + 4) = qtr_cfc (:,:,3)
273            trc2d(:,:,jp_cfc0_2d + 5) = qint_cfc(:,:,4)
274         END IF
275      END IF
276      !
277      IF( l_trdtrc ) THEN
278          DO jn = jp_cfc0, jp_cfc1
279            CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
280          END DO
281      END IF
282      !
283# if defined key_debug_medusa
284      IF(lwp) WRITE(numout,*) '   CFC - Check: nn_timing = ', nn_timing
285      CALL flush(numout)
286# endif
287      IF( nn_timing == 1 )  CALL timing_stop('trc_sms_cfc')
288      !
289   END SUBROUTINE trc_sms_cfc
290
291
292   SUBROUTINE cfc_init
293      !!---------------------------------------------------------------------
294      !!                     ***  cfc_init  *** 
295      !!
296      !! ** Purpose : sets constants for CFC model
297      !!---------------------------------------------------------------------
298      INTEGER :: jl, jn, iyear_beg, iyear_tmp
299
300      ! coefficient for CFC11
301      !----------------------
302
303      ! Solubility
304      soa(1,1) = -229.9261 
305      soa(2,1) =  319.6552
306      soa(3,1) =  119.4471
307      soa(4,1) =   -1.39165
308
309      sob(1,1) = -0.142382
310      sob(2,1) =  0.091459
311      sob(3,1) = -0.0157274
312
313      ! Schmidt number          AXY (25/04/17)
314      sca(1,1) = 3579.2       ! = 3501.8
315      sca(2,1) = -222.63      ! = -210.31
316      sca(3,1) =    7.5749    ! =    6.1851
317      sca(4,1) =   -0.14595   ! =   -0.07513
318      sca(5,1) =    0.0011874 ! = absent
319
320      ! coefficient for CFC12
321      !----------------------
322
323      ! Solubility
324      soa(1,2) = -218.0971
325      soa(2,2) =  298.9702
326      soa(3,2) =  113.8049
327      soa(4,2) =   -1.39165
328
329      sob(1,2) = -0.143566
330      sob(2,2) =  0.091015
331      sob(3,2) = -0.0153924
332
333      ! schmidt number         AXY (25/04/17)
334      sca(1,2) = 3828.1      ! = 3845.4
335      sca(2,2) = -249.86     ! = -228.95
336      sca(3,2) =    8.7603   ! =    6.1908
337      sca(4,2) =   -0.1716   ! =   -0.067430
338      sca(5,2) =    0.001408 ! = absent
339
340      ! coefficients for SF6   AXY (25/04/17)
341      !---------------------
342     
343      ! Solubility
344      soa(1,3) =  -80.0343
345      soa(2,3) =  117.232
346      soa(3,3) =   29.5817
347      soa(4,3) =    0.0
348
349      sob(1,3) =  0.0335183
350      sob(2,3) = -0.0373942
351      sob(3,3) =  0.00774862
352
353      ! Schmidt number
354      sca(1,3) = 3177.5
355      sca(2,3) = -200.57
356      sca(3,3) =    6.8865
357      sca(4,3) =   -0.13335
358      sca(5,3) =    0.0010877
359
360      !!---------------------------------------------
361      !! JPALM -- re-initialize CFC fields and diags if restart a CFC cycle,
362      !!       Or if out of P_cfc range
363      IF (simu_type==1) THEN
364         iyear_tmp = nyear - nyear_res  !! JPALM -- in our spin-up, nyear_res is 1000
365         iyear_beg = MOD( iyear_tmp , 90 )
366         !!---------------------------------------
367         IF ((iyear_beg > 77) .OR. (iyear_beg==0)) THEN
368            qtr_cfc(:,:,:) = 0._wp
369            IF(lwp) THEN
370               WRITE(numout,*) 
371               WRITE(numout,*) 'restart a CFC cycle or out of P_cfc year bounds zero --'
372               WRITE(numout,*) '                          --    set qtr_CFC = 0.00   --'
373               WRITE(numout,*) '                          --   set qint_CFC = 0.00   --'
374               WRITE(numout,*) '                          --   set trn(CFC) = 0.00   --'
375            ENDIF
376            qtr_cfc(:,:,:) = 0._wp
377            qint_cfc(:,:,:) = 0._wp
378            DO jl = 1, jp_cfc
379              jn = jp_cfc0 + jl - 1
380              trn(:,:,:,jn) = 0._wp
381              trb(:,:,:,jn) = 0._wp
382            END DO
383         ENDIF
384      !!
385      !! 2 -- Hindcast/proj
386      ELSEIF (simu_type==2) THEN
387         iyear_beg = MOD(nyear, 100)
388         IF (iyear_beg < 9)  iyear_beg = iyear_beg + 100
389         IF ((iyear_beg < 30) .OR. (iyear_beg > 107)) THEN
390            qtr_cfc(:,:,:) = 0._wp
391            IF(lwp) THEN
392               WRITE(numout,*)
393               WRITE(numout,*) 'restart a CFC cycle or out of P_cfc year bounds zero --'
394               WRITE(numout,*) '                          --    set qtr_CFC = 0.00   --'
395               WRITE(numout,*) '                          --   set qint_CFC = 0.00   --'
396               WRITE(numout,*) '                          --   set trn(CFC) = 0.00   --'
397            ENDIF
398            qtr_cfc(:,:,:) = 0._wp
399            qint_cfc(:,:,:) = 0._wp
400            DO jl = 1, jp_cfc
401              jn = jp_cfc0 + jl - 1
402              trn(:,:,:,jn) = 0._wp
403              trb(:,:,:,jn) = 0._wp
404            END DO
405         ENDIF
406      ENDIF
407
408      IF(lwp) WRITE(numout,*)
409      !
410   END SUBROUTINE cfc_init
411
412
413   INTEGER FUNCTION trc_sms_cfc_alloc()
414      !!----------------------------------------------------------------------
415      !!                     ***  ROUTINE trc_sms_cfc_alloc  ***
416      !!----------------------------------------------------------------------
417      ALLOCATE( xphem   (jpi,jpj)        ,     &
418         &      qtr_cfc (jpi,jpj,jp_cfc) ,     &
419         &      qint_cfc(jpi,jpj,jp_cfc) , STAT=trc_sms_cfc_alloc )
420         !
421      IF( trc_sms_cfc_alloc /= 0 ) CALL ctl_warn('trc_sms_cfc_alloc : failed to allocate arrays.')
422      !
423   END FUNCTION trc_sms_cfc_alloc
424
425#else
426   !!----------------------------------------------------------------------
427   !!   Dummy module                                         No CFC tracers
428   !!----------------------------------------------------------------------
429CONTAINS
430   SUBROUTINE trc_sms_cfc( kt )       ! Empty routine
431      WRITE(*,*) 'trc_sms_cfc: You should not have seen this print! error?', kt
432   END SUBROUTINE trc_sms_cfc
433#endif
434
435   !!======================================================================
436END MODULE trcsms_cfc
Note: See TracBrowser for help on using the repository browser.