MODULE p4zflx
   !!======================================================================
   !!                         ***  MODULE p4zflx  ***
   !! TOP :   PISCES CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
   !!======================================================================
   !! History :    -   !  1988-07  (E. MAIER-REIMER) Original code
   !!              -   !  1998     (O. Aumont) additions
   !!              -   !  1999     (C. Le Quere) modifications
   !!             1.0  !  2004     (O. Aumont) modifications
   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
   !!----------------------------------------------------------------------
#if defined key_pisces
   !!----------------------------------------------------------------------
   !!   'key_pisces'                                       PISCES bio-model
   !!----------------------------------------------------------------------
   !!   p4z_flx       :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
   !!----------------------------------------------------------------------
   USE oce_trc         !
   USE trp_trc
   USE sms

   IMPLICIT NONE
   PRIVATE

   PUBLIC   p4z_flx    ! called in p4zprg.F90

   !!* Substitution
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
   !! $Header:$ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE p4z_flx
      !!---------------------------------------------------------------------
      !!                     ***  ROUTINE p4z_flx  ***
      !!
      !! ** Purpose :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
      !!
      !! ** Method  : - ???
      !!---------------------------------------------------------------------
      INTEGER  ::   ji, jj, jrorr
      REAL(wp) ::   zpdtan, zttc, zws
      REAL(wp) ::   zfld, zflu, zoxy16, zflu16, zfact
      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zschmitto2, zalka, zschmittco2
      REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3
      !!---------------------------------------------------------------------

      ! -----------------------------------------------------
      !     ASSIGNATION TO EXPONENTS IN THE LISS AND MERLIVAT
      !     FORMULATION OF THE GAS EXCHANGE RATE
      ! -----------------------------------------------------

      zpdtan = raass / rdt

      ! SURFACE CHEMISTRY (PCO2 AND [H+] IN
      !     SURFACE LAYER); THE RESULT OF THIS CALCULATION
      !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2

      DO jrorr = 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  = trn(ji,jj,1,jpdic) / zfact
               zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
               zalka = trn(ji,jj,1,jptal) / zfact

               ! CALCULATE [ALK]([CO3--], [HCO3-])
               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph + 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
         END DO
      END DO


      ! --------------
      ! COMPUTE FLUXES
      ! --------------

      ! FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
      ! -------------------------------------------

      DO jj = 1, jpj
         DO ji = 1, jpi

            zttc = MIN( 35., tn(ji,jj,1) )
!!gm  optimisation & more precise computation with factorisation of the polynome
            zschmittco2 = 2073.1 - 125.62 * zttc + 3.6276 * zttc**2 - 0.043126 * zttc**3
            zws         = vatm(ji,jj)

            ! COMPUTE GAS EXCHANGE FOR CO2
            zkgco2(ji,jj) = (  0.3 * zws * zws    &
               &             + 2.5 * ( 0.5246 + zttc * ( 0.016256 + zttc * 0.00049946 ) )  )   &
               &          * SQRT( 660./ zschmittco2 )

            ! CONVERT TO m/s, and apply sea-ice cover
            zkgco2(ji,jj) = zkgco2(ji,jj) / ( 100. * 3600. )      &
# if defined key_off_degrad
               &         * facvol(ji,jj,1)      &
# endif
               &         * ( 1.- freeze(ji,jj) ) * tmask(ji,jj,1)

         END DO
      END DO

      ! COMPUTE GAS EXCHANGE COEFFICIENT FO O2 FROM Waninkhof EQUATIONS
      DO jj = 1, jpj
         DO ji = 1, jpi

          zws  = vatm(ji,jj)
          zttc = MIN( 35., tn(ji,jj,1) )
!!gm  optimisation & more precise computation with factorisation of the polynome
          zschmitto2   = 1953.4 - 128.0 * zttc + 3.9918 * zttc**2 - 0.050091 * zttc**3

          zkgo2(ji,jj) = (  0.3 * zws * zws   &
             &            + 2.5 * ( 0.5246 + zttc * ( 0.016256 + zttc * 0.00049946 ) )  )   &
             &         * SQRT( 660./ zschmitto2 )

          ! CONVERT TO m/s AND APPLY SEA ICE COVER
          zkgo2(ji,jj) = zkgo2(ji,jj) / ( 100.*3600.)        &
# if defined key_off_degrad
             &        * facvol(ji,jj,1)       &
# endif
             &        * ( 1.- freeze(ji,jj) ) *tmask(ji,jj,1)

         END DO
      END DO

      DO jj = 1, jpj
         DO ji = 1, jpi

            ! Compute CO2 flux for the sea and air
            zfld = atcco2 * tmask(ji,jj,1) * chemc(ji,jj,3) * zkgco2(ji,jj)
            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)
            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) / fse3t(ji,jj,1)

            ! Compute O2 flux 
            zoxy16 = trn(ji,jj,1,jpoxy)
            zflu16 = ( atcox * chemc(ji,jj,2) - zoxy16 ) * zkgo2(ji,jj)
            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zflu16 / fse3t(ji,jj,1)

# if defined key_trc_diaadd
            ! Save diagnostics
            trc2d(ji,jj,1) = ( zfld - zflu ) * 1000.
            trc2d(ji,jj,2) = zflu16 * 1000.
            trc2d(ji,jj,3) = zkgco2(ji,jj)
            trc2d(ji,jj,4) = atcco2 - zh2co3(ji,jj) / ( chemc(ji,jj,3) + rtrn )
# endif
         END DO
      END DO
      !
   END SUBROUTINE p4z_flx

#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 

   !!======================================================================
END MODULE  p4zflx