Index: branches/UKMO/dev_r5518_GO6M_dev/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6M_dev/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 8199)
+++ branches/UKMO/dev_r5518_GO6M_dev/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 8201)
@@ -7,4 +7,6 @@
    !!  NEMO      1.0  !  2004-03  (C. Ethe) free form + modularity
    !!            2.0  !  2007-12  (C. Ethe, G. Madec)  reorganisation
+   !!                 !  2016-06  (J. Palmieri)  update for UKESM1
+   !!                 !  2017-04  (A. Yool)  update to add SF6, fix coefficients
    !!----------------------------------------------------------------------
 #if defined key_cfc
@@ -15,4 +17,5 @@
    !!   cfc_init     :  sets constants for CFC surface forcing computation
    !!----------------------------------------------------------------------
+   USE dom_oce       ! ocean space and time domain
    USE oce_trc       ! Ocean variables
    USE par_trc       ! TOP parameters
@@ -31,4 +34,6 @@
    INTEGER , PUBLIC            ::   jpyear         ! Number of years read in CFC1112 file
    INTEGER , PUBLIC            ::   ndate_beg      ! initial calendar date (aammjj) for CFC
+   INTEGER , PUBLIC            ::   simu_type      ! Kind of simulation: 1- Spin-up 
+                                                   !                     2- Hindcast/projection
    INTEGER , PUBLIC            ::   nyear_res      ! restoring time constant (year)
    INTEGER , PUBLIC            ::   nyear_beg      ! initial year (aa) 
@@ -40,7 +45,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   patm     ! atmospheric function
 
-   REAL(wp), DIMENSION(4,2) ::   soa   ! coefficient for solubility of CFC [mol/l/atm]
-   REAL(wp), DIMENSION(3,2) ::   sob   !    "               "
-   REAL(wp), DIMENSION(4,2) ::   sca   ! coefficients for schmidt number in degre Celcius
+   REAL(wp), DIMENSION(4,3) ::   soa   ! coefficient for solubility of CFC [mol/l/atm]
+   REAL(wp), DIMENSION(3,3) ::   sob   !    "               "
+   REAL(wp), DIMENSION(5,3) ::   sca   ! coefficients for schmidt number in degre Celcius
       
    !                          ! coefficients for conversion
@@ -79,8 +84,8 @@
       !
       INTEGER  ::   ji, jj, jn, jl, jm, js
-      INTEGER  ::   iyear_beg, iyear_end
+      INTEGER  ::   iyear_beg, iyear_end, iyear_tmp
       INTEGER  ::   im1, im2, ierr
       REAL(wp) ::   ztap, zdtap        
-      REAL(wp) ::   zt1, zt2, zt3, zv2
+      REAL(wp) ::   zt1, zt2, zt3, zt4, zv2
       REAL(wp) ::   zsol      ! solubility
       REAL(wp) ::   zsch      ! schmidt number 
@@ -103,5 +108,38 @@
       ! Temporal interpolation
       ! ----------------------
-      iyear_beg = nyear - 1900
+      !! JPALM -- 15-06-2016 -- define 2 kinds of CFC run:
+      !!                     1- the SPIN-UP and 2- Hindcast/Projections
+      !!                     -- main difference is the way to define the year of
+      !!                     simulation, that determine the atm pCFC.
+      !!                     1-- Spin-up: our atm forcing is of 30y we cycle on.
+      !!                     So we do 90y CFC cycles to be in good
+      !!                     correspondence with the atmosphere
+      !!                     2-- Hindcast/proj, instead of nyear-1900 we keep
+      !!                     the 2 last digit, and enable 3 cycle from 1800 to 2100.  
+      !!----------------------------------------------------------------------
+      IF (simu_type==1) THEN
+         !! 1 -- SPIN-UP
+         iyear_tmp = nyear - nyear_res  !! JPALM -- in our spin-up, nyear_res is 1000
+         iyear_beg = MOD( iyear_tmp , 90 )
+         !! JPALM -- the pCFC file only got 78 years.
+         !!       So if iyear_beg > 78 then we set pCFC to 0
+         !!             iyear_beg = 0 as well -- must try to avoid obvious problems
+         !!             as Pcfc is set to 0.00 up to year 32, let set iyear_beg to year 10
+         !!          else, must add 30 to iyear_beg to match with P_cfc indices
+         !!---------------------------------------
+         IF ((iyear_beg > 77) .OR. (iyear_beg==0)) THEN
+            iyear_beg = 10
+         ELSE 
+            iyear_beg = iyear_beg + 30
+         ENDIF
+      ELSEIF (simu_type==2) THEN
+         !! 2 -- Hindcast/proj
+         iyear_beg = MOD(nyear, 100)
+         IF (iyear_beg < 20)  iyear_beg = iyear_beg + 100
+         !! JPALM -- Same than previously, if iyear_beg is out of P_cfc range,
+         !!       we want to set p_CFC to 0.00 --> set iyear_beg = 10
+         IF ((iyear_beg < 30) .OR. (iyear_beg > 115)) iyear_beg = 10             
+      ENDIF
+      !!
       IF ( nmonth <= 6 ) THEN
          iyear_beg = iyear_beg - 1
@@ -152,10 +190,13 @@
                zt2  = zt1 * zt1 
                zt3  = zt1 * zt2
-               zsch = sca(1,jl) + sca(2,jl) * zt1 + sca(3,jl) * zt2 + sca(4,jl) * zt3
+               zt4  = zt1 * zt3
+               zsch = sca(1,jl) + sca(2,jl) * zt1 + sca(3,jl) * zt2 + sca(4,jl) * zt3 + sca(5,jl) * zt4
 
                !    speed transfert : formulae of wanninkhof 1992
                zv2     = wndm(ji,jj) * wndm(ji,jj)
                zsch    = zsch / 660.
-               zak_cfc = ( 0.39 * xconv2 * zv2 / SQRT(zsch) ) * tmask(ji,jj,1)
+               ! AXY (25/04/17): OMIP protocol specifies lower Wanninkhof (2014) value
+               ! zak_cfc = ( 0.39 * xconv2 * zv2 / SQRT(zsch) ) * tmask(ji,jj,1)
+               zak_cfc = ( 0.251 * xconv2 * zv2 / SQRT(zsch) ) * tmask(ji,jj,1)
 
                ! Input function  : speed *( conc. at equil - concen at surface )
@@ -176,22 +217,59 @@
          !                                                  !----------------!
       END DO                                                !  end CFC loop  !
-      !
-      IF( lrst_trc ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'trc_sms_cfc : cumulated input function fields written in ocean restart file ',   &
-            &                    'at it= ', kt,' date= ', ndastp
-         IF(lwp) WRITE(numout,*) '~~~~'
-         DO jn = jp_cfc0, jp_cfc1
-            CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
-         END DO
-      ENDIF                                            
+         !
+      IF( kt == nittrc000 ) THEN
+         DO jl = 1, jp_cfc   
+             WRITE(NUMOUT,*) ' '
+             WRITE(NUMOUT,*) 'CFC interpolation verification '  !! Jpalm  
+             WRITE(NUMOUT,*) '################################## '
+             WRITE(NUMOUT,*) ' '
+               if (jl.EQ.1) then
+                   WRITE(NUMOUT,*) 'Traceur = CFC11: '
+               elseif (jl.EQ.2) then
+                   WRITE(NUMOUT,*) 'Traceur = CFC12: '
+               elseif (jl.EQ.3) then
+                   WRITE(NUMOUT,*) 'Traceur = SF6: '
+               endif
+             WRITE(NUMOUT,*) 'nyear    = ', nyear
+             WRITE(NUMOUT,*) 'nmonth   = ', nmonth
+             WRITE(NUMOUT,*) 'iyear_beg= ', iyear_beg
+             WRITE(NUMOUT,*) 'iyear_end= ', iyear_end
+             WRITE(NUMOUT,*) 'p_cfc(iyear_beg)= ',p_cfc(iyear_beg, 1, jl)
+             WRITE(NUMOUT,*) 'p_cfc(iyear_end)= ',p_cfc(iyear_end, 1, jl)
+             WRITE(NUMOUT,*) 'Im1= ',im1
+             WRITE(NUMOUT,*) 'Im2= ',im2
+             WRITE(NUMOUT,*) 'zpp_cfc = ',zpp_cfc
+             WRITE(NUMOUT,*) ' '
+         END DO  
+# if defined key_debug_medusa
+         CALL flush(numout)
+# endif
+      ENDIF
+        !
+      !IF( lrst_trc ) THEN
+      !   IF(lwp) WRITE(numout,*)
+      !   IF(lwp) WRITE(numout,*) 'trc_sms_cfc : cumulated input function fields written in ocean restart file ',   &
+      !      &                    'at it= ', kt,' date= ', ndastp
+      !   IF(lwp) WRITE(numout,*) '~~~~'
+      !   DO jn = jp_cfc0, jp_cfc1
+      !      CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
+      !   END DO
+      !ENDIF                                            
       !
       IF( lk_iomput ) THEN
-         CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
-         CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
+         IF  (iom_use("qtrCFC11"))  CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
+         IF  (iom_use("qintCFC11")) CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
+         IF  (iom_use("qtrCFC12"))  CALL iom_put( "qtrCFC12"  , qtr_cfc (:,:,2) )
+         IF  (iom_use("qintCFC12")) CALL iom_put( "qintCFC12" , qint_cfc(:,:,2) )
+         IF  (iom_use("qtrSF6"))    CALL iom_put( "qtrSF6"    , qtr_cfc (:,:,3) )
+         IF  (iom_use("qintSF6"))   CALL iom_put( "qintSF6"   , qint_cfc(:,:,3) )
       ELSE
          IF( ln_diatrc ) THEN
             trc2d(:,:,jp_cfc0_2d    ) = qtr_cfc (:,:,1)
             trc2d(:,:,jp_cfc0_2d + 1) = qint_cfc(:,:,1)
+            trc2d(:,:,jp_cfc0_2d + 2) = qtr_cfc (:,:,2)
+            trc2d(:,:,jp_cfc0_2d + 3) = qint_cfc(:,:,2)
+            trc2d(:,:,jp_cfc0_2d + 4) = qtr_cfc (:,:,3)
+            trc2d(:,:,jp_cfc0_2d + 5) = qint_cfc(:,:,3)
          END IF
       END IF
@@ -203,4 +281,8 @@
       END IF
       !
+# if defined key_debug_medusa
+      IF(lwp) WRITE(numout,*) '   CFC - Check: nn_timing = ', nn_timing
+      CALL flush(numout)
+# endif
       IF( nn_timing == 1 )  CALL timing_stop('trc_sms_cfc')
       !
@@ -214,5 +296,5 @@
       !! ** Purpose : sets constants for CFC model
       !!---------------------------------------------------------------------
-      INTEGER :: jn
+      INTEGER :: jl, jn, iyear_beg, iyear_tmp
 
       ! coefficient for CFC11 
@@ -223,15 +305,16 @@
       soa(2,1) =  319.6552
       soa(3,1) =  119.4471
-      soa(4,1) =  -1.39165
-
-      sob(1,1) =  -0.142382
-      sob(2,1) =   0.091459
-      sob(3,1) =  -0.0157274
-
-      ! Schmidt number 
-      sca(1,1) = 3501.8
-      sca(2,1) = -210.31
-      sca(3,1) =  6.1851
-      sca(4,1) = -0.07513
+      soa(4,1) =   -1.39165
+
+      sob(1,1) = -0.142382
+      sob(2,1) =  0.091459
+      sob(3,1) = -0.0157274
+
+      ! Schmidt number          AXY (25/04/17)
+      sca(1,1) = 3579.2       ! = 3501.8
+      sca(2,1) = -222.63      ! = -210.31
+      sca(3,1) =    7.5749    ! =    6.1851
+      sca(4,1) =   -0.14595   ! =   -0.07513
+      sca(5,1) =    0.0011874 ! = absent
 
       ! coefficient for CFC12 
@@ -242,25 +325,79 @@
       soa(2,2) =  298.9702
       soa(3,2) =  113.8049
-      soa(4,2) =  -1.39165
-
-      sob(1,2) =  -0.143566
-      sob(2,2) =   0.091015
-      sob(3,2) =  -0.0153924
-
-      ! schmidt number 
-      sca(1,2) =  3845.4 
-      sca(2,2) =  -228.95
-      sca(3,2) =  6.1908 
-      sca(4,2) =  -0.067430
-
-      IF( ln_rsttr ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) ' Read specific variables from CFC model '
-         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
-         !
-         DO jn = jp_cfc0, jp_cfc1
-            CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) ) 
-         END DO
+      soa(4,2) =   -1.39165
+
+      sob(1,2) = -0.143566
+      sob(2,2) =  0.091015
+      sob(3,2) = -0.0153924
+
+      ! schmidt number         AXY (25/04/17)
+      sca(1,2) = 3828.1      ! = 3845.4 
+      sca(2,2) = -249.86     ! = -228.95
+      sca(3,2) =    8.7603   ! =    6.1908 
+      sca(4,2) =   -0.1716   ! =   -0.067430
+      sca(5,2) =    0.001408 ! = absent
+
+      ! coefficients for SF6   AXY (25/04/17)
+      !---------------------
+      
+      ! Solubility
+      soa(1,3) =  -80.0343
+      soa(2,3) =  117.232
+      soa(3,3) =   29.5817
+      soa(4,3) =    0.0
+
+      sob(1,3) =  0.0335183
+      sob(2,3) = -0.0373942
+      sob(3,3) =  0.00774862
+
+      ! Schmidt number
+      sca(1,3) = 3177.5
+      sca(2,3) = -200.57
+      sca(3,3) =    6.8865
+      sca(4,3) =   -0.13335
+      sca(5,3) =    0.0010877
+
+      !!---------------------------------------------
+      !! JPALM -- re-initialize CFC fields and diags if restart a CFC cycle,
+      !!       Or if out of P_cfc range
+      IF (simu_type==1) THEN
+         iyear_tmp = nyear - nyear_res  !! JPALM -- in our spin-up, nyear_res is 1000
+         iyear_beg = MOD( iyear_tmp , 90 )
+         !!---------------------------------------
+         IF ((iyear_beg > 77) .OR. (iyear_beg==0)) THEN
+            qtr_cfc(:,:,:) = 0._wp
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'restart a CFC cycle or out of P_cfc year bounds zero --'
+               WRITE(numout,*) '                          --    set qtr_CFC = 0.00   --'
+               WRITE(numout,*) '                          --   set qint_CFC = 0.00   --'
+               WRITE(numout,*) '                          --   set trn(CFC) = 0.00   --'
+            ENDIF
+            qtr_cfc(:,:,:) = 0._wp
+            qint_cfc(:,:,:) = 0._wp
+            trn(:,:,:,jp_cfc0:jp_cfc1) = 0._wp
+            trb(:,:,:,jp_cfc0:jp_cfc1) = 0._wp
+         ENDIF
+      !!
+      !! 2 -- Hindcast/proj
+      ELSEIF (simu_type==2) THEN
+         iyear_beg = MOD(nyear, 100)
+         IF (iyear_beg < 20)  iyear_beg = iyear_beg + 100
+         IF ((iyear_beg < 30) .OR. (iyear_beg > 115)) THEN
+            qtr_cfc(:,:,:) = 0._wp
+            IF(lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) 'restart a CFC cycle or out of P_cfc year bounds zero --'
+               WRITE(numout,*) '                          --    set qtr_CFC = 0.00   --'
+               WRITE(numout,*) '                          --   set qint_CFC = 0.00   --'
+               WRITE(numout,*) '                          --   set trn(CFC) = 0.00   --'
+            ENDIF
+            qtr_cfc(:,:,:) = 0._wp
+            qint_cfc(:,:,:) = 0._wp
+            trn(:,:,:,jp_cfc0:jp_cfc1) = 0._wp
+            trb(:,:,:,jp_cfc0:jp_cfc1) = 0._wp
+         ENDIF
       ENDIF
+
       IF(lwp) WRITE(numout,*)
       !
