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

Last change on this file since 12240 was 12240, checked in by frrh, 11 months ago

Manual change applied to TOP_SRC/CFC/trcsms_cfc.F90 to reflect
change from branches/UKMO/dev_r5518_GO6_package_text_diagnostics
revision 12228.

This is a result of an unexplained failure of the merge command for this
routine only, while the two other files involved, icbclv.F90 and
stpctl.F90 worked as expected.

The symptoms indicated that the change had been applied locally in the working
copy of dev_r5518_GO6_package, but svn refused to acknowledge that
any change had been made and failed to apply the change at commit time.

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