Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 7381)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 7391)
@@ -11,8 +11,4 @@
    !!                  !  2011-02  (J. Simeon, J. Orr) Include total atm P correction 
    !!----------------------------------------------------------------------
-#if defined key_pisces
-   !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
-   !!----------------------------------------------------------------------
    !!   p4z_flx       :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
    !!   p4z_flx_init  :   Read the namelist
@@ -26,7 +22,4 @@
    USE iom                          !  I/O manager
    USE fldread                      !  read input fields
-#if defined key_cpl_carbon_cycle
-   USE sbc_oce, ONLY :  atm_co2     !  atmospheric pCO2               
-#endif
 
    IMPLICIT NONE
@@ -48,11 +41,11 @@
 
    !                               !!* nampisatm namelist (Atmospheric PRessure) *
-   LOGICAL, PUBLIC ::   ln_presatm  !: ref. pressure: global mean Patm (F) or a constant (F)
-
-   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:)  ::  patm      ! atmospheric pressure at kt                 [N/m2]
-   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)    ::  sf_patm   ! structure of input fields (file informations, fields read)
-
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: oce_co2   !: ocean carbon flux 
+   LOGICAL, PUBLIC ::   ln_presatm     !: ref. pressure: global mean Patm (F) or a constant (F)
+   LOGICAL, PUBLIC ::   ln_presatmco2  !: accounting for spatial atm CO2 in the compuation of carbon flux (T) or not (F)
+
+   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) ::  patm      ! atmospheric pressure at kt                 [N/m2]
+   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)   ::  sf_patm   ! structure of input fields (file informations, fields read)
+   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)   ::  sf_atmco2 ! structure of input fields (file informations, fields read)
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: satmco2   !: atmospheric pco2 
 
@@ -74,5 +67,4 @@
       !! ** Method  : 
       !!              - Include total atm P correction via Esbensen & Kushnir (1981) 
-      !!              - Pressure correction NOT done for key_cpl_carbon_cycle
       !!              - Remove Wanninkhof chemical enhancement;
       !!              - Add option for time-interpolation of atcco2.txt  
@@ -85,5 +77,5 @@
       REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
       REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff
-      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2
+      REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2
       REAL(wp) ::   zyr_dec, zdco2dt
       CHARACTER (len=25) :: charout
@@ -100,7 +92,7 @@
       !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
 
-      IF( kt /= nit000 .AND. knt == 1 ) CALL p4z_patm( kt )    ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
-
-      IF( ln_co2int ) THEN 
+      IF( kt /= nit000 .AND. .NOT.l_co2cpl .AND. knt == 1 ) CALL p4z_patm( kt )    ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
+
+      IF( ln_co2int .AND. .NOT.ln_presatmco2 .AND. .NOT.l_co2cpl ) THEN 
          ! Linear temporal interpolation  of atmospheric pco2.  atcco2.txt has annual values.
          ! Caveats: First column of .txt must be in years, decimal  years preferably. 
@@ -116,33 +108,16 @@
       ENDIF
 
-#if defined key_cpl_carbon_cycle
-      satmco2(:,:) = atm_co2(:,:)
-#endif
-
-      DO jm = 1, 10
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
-               zbot  = borat(ji,jj,1)
-               zfact = rhop(ji,jj,1) / 1000. + rtrn
-               zdic  = trb(ji,jj,1,jpdic) / zfact
-               zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
-               zalka = trb(ji,jj,1,jptal) / zfact
-
-               ! CALCULATE [ALK]([CO3--], [HCO3-])
-               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph / aphscale(ji,jj,1)    &
-               &       + zbot / ( 1.+ zph / akb3(ji,jj,1) )  )
-
-               ! CALCULATE [H+] AND [H2CO3]
-               zah2   = SQRT(  (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1)   &
-                  &                                        / ak13(ji,jj,1) ) * ( 2.* zdic - zalk )  )
-               zah2   = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 )
-               zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact
-               hi(ji,jj,1)   = zah2 * zfact
-            END DO
+      IF( l_co2cpl )   satmco2(:,:) = atm_co2(:,:)
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
+            zfact = rhop(ji,jj,1) / 1000. + rtrn
+            zdic  = trb(ji,jj,1,jpdic)
+            zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
+            ! CALCULATE [H2CO3]
+            zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
          END DO
       END DO
-
 
       ! --------------
@@ -167,7 +142,4 @@
             zkgwan = 0.251 * zws
             zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
-# if defined key_degrad
-            zkgwan = zkgwan * facvol(ji,jj,1)
-#endif 
             ! compute gas exchange for CO2 and O2
             zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 )
@@ -176,8 +148,9 @@
       END DO
 
+
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ztkel     = tsn(ji,jj,1,jp_tem) + 273.15
-            zsal      = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35.
+            ztkel = tempis(ji,jj,1) + 273.15
+            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
             zvapsw    = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
             zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
@@ -232,5 +205,5 @@
          ENDIF
          IF( iom_use( "Dpo2" ) )  THEN
-           zw2d(:,:) = ( atcox * patm(:,:) - atcox * trn(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1)
+           zw2d(:,:) = ( atcox * patm(:,:) - atcox * trb(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1)
            CALL iom_put( "Dpo2"  , zw2d )
          ENDIF
@@ -239,11 +212,4 @@
          !
          CALL wrk_dealloc( jpi, jpj, zw2d )
-      ELSE
-         IF( ln_diatrc ) THEN
-            trc2d(:,:,jp_pcs0_2d    ) = oce_co2(:,:) / e1e2t(:,:) * rfact2r 
-            trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 
-            trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) 
-            trc2d(:,:,jp_pcs0_2d + 3) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
-         ENDIF
       ENDIF
       !
@@ -287,5 +253,5 @@
          WRITE(numout,*) ' '
       ENDIF
-      IF( .NOT.ln_co2int ) THEN
+     IF( .NOT.ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
          IF(lwp) THEN                         ! control print
             WRITE(numout,*) '    Constant Atmospheric pCO2 value  atcco2    =', atcco2
@@ -293,5 +259,5 @@
          ENDIF
          satmco2(:,:)  = atcco2      ! Initialisation of atmospheric pco2
-      ELSE
+      ELSEIF( ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
          IF(lwp)  THEN
             WRITE(numout,*) '    Atmospheric pCO2 value  from file clname      =', TRIM( clname )
@@ -315,5 +281,16 @@
          END DO
          CLOSE(numco2)
-      ENDIF
+      ELSEIF( .NOT.ln_co2int .AND. ln_presatmco2 ) THEN
+         IF(lwp)  THEN
+            WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
+            WRITE(numout,*) ' '
+         ENDIF
+      ELSE
+         IF(lwp)  THEN
+            WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
+            WRITE(numout,*) ' '
+         ENDIF
+      ENDIF
+
       !
       oce_co2(:,:)  = 0._wp                ! Initialization of Flux of Carbon
@@ -341,6 +318,7 @@
       CHARACTER(len=100) ::  cn_dir   ! Root directory for location of ssr files
       TYPE(FLD_N)        ::  sn_patm  ! informations about the fields to be read
-      !!
-      NAMELIST/nampisatm/ ln_presatm, sn_patm, cn_dir
+      TYPE(FLD_N)        ::  sn_atmco2 ! informations about the fields to be read
+      !!
+      NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir
 
       !                                         ! ----------------------- !
@@ -361,4 +339,5 @@
             WRITE(numout,*) '   Namelist nampisatm : Atmospheric Pressure as external forcing'
             WRITE(numout,*) '      constant atmopsheric pressure (F) or from a file (T)  ln_presatm = ', ln_presatm
+            WRITE(numout,*) '      spatial atmopsheric CO2 for flux calcs  ln_presatmco2 = ', ln_presatmco2
             WRITE(numout,*)
          ENDIF
@@ -373,4 +352,13 @@
          ENDIF
          !                                         
+         IF( ln_presatmco2 ) THEN
+            ALLOCATE( sf_atmco2(1), STAT=ierr )           !* allocate and fill sf_atmco2 (forcing structure) with sn_atmco2
+            IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_atmco2 structure' )
+            !
+            CALL fld_fill( sf_atmco2, (/ sn_atmco2 /), cn_dir, 'p4z_flx', 'Atmospheric co2 partial pressure ', 'nampisatm' )
+                                   ALLOCATE( sf_atmco2(1)%fnow(jpi,jpj,1)   )
+            IF( sn_atmco2%ln_tint )  ALLOCATE( sf_atmco2(1)%fdta(jpi,jpj,1,2) )
+         ENDIF
+         !
          IF( .NOT.ln_presatm )   patm(:,:) = 1.e0    ! Initialize patm if no reading from a file
          !
@@ -382,26 +370,23 @@
       ENDIF
       !
+      IF( ln_presatmco2 ) THEN
+         CALL fld_read( kt, 1, sf_atmco2 )               !* input atmco2 provided at kt + 1/2
+         satmco2(:,:) = sf_atmco2(1)%fnow(:,:,1)                        ! atmospheric pressure
+      ELSE
+         satmco2(:,:) = atcco2    ! Initialize atmco2 if no reading from a file
+      ENDIF
+      !
    END SUBROUTINE p4z_patm
 
+
    INTEGER FUNCTION p4z_flx_alloc()
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_flx_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( oce_co2(jpi,jpj), satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc )
+      ALLOCATE( satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc )
       !
       IF( p4z_flx_alloc /= 0 )   CALL ctl_warn('p4z_flx_alloc : failed to allocate arrays')
       !
    END FUNCTION p4z_flx_alloc
-
-#else
-   !!======================================================================
-   !!  Dummy module :                                   No PISCES bio-model
-   !!======================================================================
-CONTAINS
-   SUBROUTINE p4z_flx( kt )                   ! Empty routine
-      INTEGER, INTENT( in ) ::   kt
-      WRITE(*,*) 'p4z_flx: You should not have seen this print! error?', kt
-   END SUBROUTINE p4z_flx
-#endif 
 
    !!======================================================================
