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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/CFC – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90 @ 8442

Last change on this file since 8442 was 8442, checked in by frrh, 7 years ago

Commit changes relating to Met Office GMED ticket 340 for the
tidying of MEDUSA related code and debugging statements in the TOP code.

Only code introduced at revision 8434 of branch
http://fcm3/projects/NEMO.xm/log/branches/NERC/dev_r5518_GO6_split_trcbiomedusa
is included here, all previous revisions of that branch having been dealt with
under GMED ticket 339.

File size: 18.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   !!                 !  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 < 20)  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 > 115)) 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  (iom_use("qtrCFC11"))  CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
260      IF  (iom_use("qintCFC11")) CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
261      IF  (iom_use("qtrCFC12"))  CALL iom_put( "qtrCFC12"  , qtr_cfc (:,:,2) )
262      IF  (iom_use("qintCFC12")) CALL iom_put( "qintCFC12" , qint_cfc(:,:,2) )
263      IF  (iom_use("qtrSF6"))    CALL iom_put( "qtrSF6"    , qtr_cfc (:,:,3) )
264      IF  (iom_use("qintSF6"))   CALL iom_put( "qintSF6"   , qint_cfc(:,:,3) )
265      !
266      IF( l_trdtrc ) THEN
267          DO jn = jp_cfc0, jp_cfc1
268            CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
269          END DO
270      END IF
271      !
272# if defined key_debug_medusa
273      IF(lwp) WRITE(numout,*) '   CFC - Check: nn_timing = ', nn_timing
274      CALL flush(numout)
275# endif
276      IF( nn_timing == 1 )  CALL timing_stop('trc_sms_cfc')
277      !
278   END SUBROUTINE trc_sms_cfc
279
280
281   SUBROUTINE cfc_init
282      !!---------------------------------------------------------------------
283      !!                     ***  cfc_init  *** 
284      !!
285      !! ** Purpose : sets constants for CFC model
286      !!---------------------------------------------------------------------
287      INTEGER :: jl, jn, iyear_beg, iyear_tmp
288
289      ! coefficient for CFC11
290      !----------------------
291
292      ! Solubility
293      soa(1,1) = -229.9261 
294      soa(2,1) =  319.6552
295      soa(3,1) =  119.4471
296      soa(4,1) =   -1.39165
297
298      sob(1,1) = -0.142382
299      sob(2,1) =  0.091459
300      sob(3,1) = -0.0157274
301
302      ! Schmidt number          AXY (25/04/17)
303      sca(1,1) = 3579.2       ! = 3501.8
304      sca(2,1) = -222.63      ! = -210.31
305      sca(3,1) =    7.5749    ! =    6.1851
306      sca(4,1) =   -0.14595   ! =   -0.07513
307      sca(5,1) =    0.0011874 ! = absent
308
309      ! coefficient for CFC12
310      !----------------------
311
312      ! Solubility
313      soa(1,2) = -218.0971
314      soa(2,2) =  298.9702
315      soa(3,2) =  113.8049
316      soa(4,2) =   -1.39165
317
318      sob(1,2) = -0.143566
319      sob(2,2) =  0.091015
320      sob(3,2) = -0.0153924
321
322      ! schmidt number         AXY (25/04/17)
323      sca(1,2) = 3828.1      ! = 3845.4
324      sca(2,2) = -249.86     ! = -228.95
325      sca(3,2) =    8.7603   ! =    6.1908
326      sca(4,2) =   -0.1716   ! =   -0.067430
327      sca(5,2) =    0.001408 ! = absent
328
329      ! coefficients for SF6   AXY (25/04/17)
330      !---------------------
331     
332      ! Solubility
333      soa(1,3) =  -80.0343
334      soa(2,3) =  117.232
335      soa(3,3) =   29.5817
336      soa(4,3) =    0.0
337
338      sob(1,3) =  0.0335183
339      sob(2,3) = -0.0373942
340      sob(3,3) =  0.00774862
341
342      ! Schmidt number
343      sca(1,3) = 3177.5
344      sca(2,3) = -200.57
345      sca(3,3) =    6.8865
346      sca(4,3) =   -0.13335
347      sca(5,3) =    0.0010877
348
349      !!---------------------------------------------
350      !! JPALM -- re-initialize CFC fields and diags if restart a CFC cycle,
351      !!       Or if out of P_cfc range
352      IF (simu_type==1) THEN
353         iyear_tmp = nyear - nyear_res  !! JPALM -- in our spin-up, nyear_res is 1000
354         iyear_beg = MOD( iyear_tmp , 90 )
355         !!---------------------------------------
356         IF ((iyear_beg > 77) .OR. (iyear_beg==0)) THEN
357            qtr_cfc(:,:,:) = 0._wp
358            IF(lwp) THEN
359               WRITE(numout,*) 
360               WRITE(numout,*) 'restart a CFC cycle or out of P_cfc year bounds zero --'
361               WRITE(numout,*) '                          --    set qtr_CFC = 0.00   --'
362               WRITE(numout,*) '                          --   set qint_CFC = 0.00   --'
363               WRITE(numout,*) '                          --   set trn(CFC) = 0.00   --'
364            ENDIF
365            qtr_cfc(:,:,:) = 0._wp
366            qint_cfc(:,:,:) = 0._wp
367            trn(:,:,:,jp_cfc0:jp_cfc1) = 0._wp
368            trb(:,:,:,jp_cfc0:jp_cfc1) = 0._wp
369         ENDIF
370      !!
371      !! 2 -- Hindcast/proj
372      ELSEIF (simu_type==2) THEN
373         iyear_beg = MOD(nyear, 100)
374         IF (iyear_beg < 20)  iyear_beg = iyear_beg + 100
375         IF ((iyear_beg < 30) .OR. (iyear_beg > 115)) THEN
376            qtr_cfc(:,:,:) = 0._wp
377            IF(lwp) THEN
378               WRITE(numout,*)
379               WRITE(numout,*) 'restart a CFC cycle or out of P_cfc year bounds zero --'
380               WRITE(numout,*) '                          --    set qtr_CFC = 0.00   --'
381               WRITE(numout,*) '                          --   set qint_CFC = 0.00   --'
382               WRITE(numout,*) '                          --   set trn(CFC) = 0.00   --'
383            ENDIF
384            qtr_cfc(:,:,:) = 0._wp
385            qint_cfc(:,:,:) = 0._wp
386            trn(:,:,:,jp_cfc0:jp_cfc1) = 0._wp
387            trb(:,:,:,jp_cfc0:jp_cfc1) = 0._wp
388         ENDIF
389      ENDIF
390
391      IF(lwp) WRITE(numout,*)
392      !
393   END SUBROUTINE cfc_init
394
395
396   INTEGER FUNCTION trc_sms_cfc_alloc()
397      !!----------------------------------------------------------------------
398      !!                     ***  ROUTINE trc_sms_cfc_alloc  ***
399      !!----------------------------------------------------------------------
400      ALLOCATE( xphem   (jpi,jpj)        ,     &
401         &      qtr_cfc (jpi,jpj,jp_cfc) ,     &
402         &      qint_cfc(jpi,jpj,jp_cfc) , STAT=trc_sms_cfc_alloc )
403         !
404      IF( trc_sms_cfc_alloc /= 0 ) CALL ctl_warn('trc_sms_cfc_alloc : failed to allocate arrays.')
405      !
406   END FUNCTION trc_sms_cfc_alloc
407
408#else
409   !!----------------------------------------------------------------------
410   !!   Dummy module                                         No CFC tracers
411   !!----------------------------------------------------------------------
412CONTAINS
413   SUBROUTINE trc_sms_cfc( kt )       ! Empty routine
414      WRITE(*,*) 'trc_sms_cfc: You should not have seen this print! error?', kt
415   END SUBROUTINE trc_sms_cfc
416#endif
417
418   !!======================================================================
419END MODULE trcsms_cfc
Note: See TracBrowser for help on using the repository browser.