source: trunk/NEMO/TOP_SRC/SMS/p4zflx.F @ 336

CCC$Header$
CCC  TOP 1.0 , LOCEAN-IPSL (2005)
C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
C ---------------------------------------------------------------------------
CDIR$ LIST
      SUBROUTINE p4zflx
#if defined key_passivetrc && defined key_trc_pisces
CCC---------------------------------------------------------------------
CCC
CCC          ROUTINE p4zflx : PISCES MODEL
CCC          *****************************
CCC
CCC
CC     PURPOSE.
CC     --------
CC          *P4ZFLX* CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
CC
CC     EXTERNALS.
CC     ----------
CC          NONE.
CC
CC   MODIFICATIONS:
CC   --------------
CC      original      : 1988-07 E. MAIER-REIMER      MPI HAMBURG
CC      additions     : 1998    O. Aumont
CC      modifications : 1999    C. Le Quere
CC      modifications : 2004    O. Aumont
CC     -----------------------------------------------------------------
CC  parameters and commons
CC ======================
CDIR$ NOLIST
      USE oce_trc
      USE trp_trc
      USE sms
      IMPLICIT NONE
CDIR$ LIST
CC----------------------------------------------------------------------
CC local declarations
CC ==================
C
      INTEGER ji, jj, krorr
      REAL zexp1, zexp2
      REAL a1, a2, a3, b2, b3, ttc, ws, alpco2
      REAL fld, flu, oxy16, flu16
      REAL zph,ah2,zbot,zdic,zalk,schmitt, zrhocd
      REAL zwind(jpi,jpj)

CC
CC----------------------------------------------------------------------
CC statement functions
CC ===================
CDIR$ NOLIST
#include "domzgr_substitute.h90"
CDIR$ LIST
C
C
C  1. ASSIGNATION TO EXPONENTS IN THE LISS AND MERLIVAT
C     FORMULATION OF THE GAS EXCHANGE RATE
c -----------------------------------------------------
C

      zexp1 = -2./3.
      zexp2 = -1./2.  
      a1    = 0.17
      a2    = 2.85
      a3    = 5.90
      b2    = 9.65
      b3    = 49.3  
 
      zrhocd = 1.3*1.3e-3        
      DO jj = 1, jpj
         DO ji = 1 , jpi
            IF (igaswind .EQ. 0 ) then
               zwind(ji,jj) = sqrt(taux(ji,jj)**2+tauy(ji,jj)**2)
     $                            /zrhocd
            ELSE
               zwind(ji,jj) = vatm(ji,jj)
            ENDIF 
         END DO
      END DO
C
C* 1.1 SURFACE CHEMISTRY (PCO2 AND [H+] IN
C     SURFACE LAYER); THE RESULT OF THIS CALCULATION
C     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
C ---------------------------------------------------
C
      DO krorr = 1,10
C
        DO jj = 1,jpj
          DO ji = 1,jpi
C
C* 1.2 DUMMY VARIABLES FOR DIC, H+, AND BORATE
C --------------------------------------------
C
            zbot = borat(ji,jj,1)
            zdic  = trn(ji,jj,1,jpdic)
            zph = max(hi(ji,jj,1),1.E-10)
C
C* 1.3 CALCULATE [ALK]([CO3--], [HCO3-])
C ------------------------------------
C
            zalk=trn(ji,jj,1,jptal)-
     &          (akw3(ji,jj,1)/zph-zph+zbot/(1.+zph/akb3(ji,jj,1)))
C
C* 1.4 CALCULATE [H+] AND [H2CO3]
C -----------------------------------------
C
            ah2=sqrt((zdic-zalk)**2+4*(zalk*ak23(ji,jj,1)
     &          /ak13(ji,jj,1))*(2*zdic-zalk))
            ah2=0.5*ak13(ji,jj,1)/zalk*((zdic-zalk)+ah2)
            h2co3(ji,jj) = (2*zdic-zalk)/(2.+ak13(ji,jj,1)/ah2)
            hi(ji,jj,1)  = ah2
          END DO
        END DO
      END DO
C
C
C 2. COMPUTE FLUXES
C --------------
C
C 2.1 FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
C -------------------------------------------
C
      DO jj = 1,jpj
        DO ji = 1,jpi

           ws = zwind(ji,jj)
          ttc = min(35.,tn(ji,jj,1))
          schmitt= 2073.1-125.62*ttc+3.6276*ttc**2-0.043126*ttc**3
C
C 2.2 COMPUTE GAS EXCHANGE FOR CO2
C --------------------------------
C
          kgwanin(ji,jj) = (0.3*ws*ws + 2.5*(0.5246+ttc*(0.016256+
     &       ttc*0.00049946)))*sqrt(660./schmitt)
C
C 2.3 CONVERT TO M/S
C ------------------
C
          kgwanin(ji,jj) = kgwanin(ji,jj)/100./3600.
C
C 2.4 convert to mol/m2/s/uatm, alpco2(chemc(ji,jj,1)) is in 
C      mol/L/uatm and apply ice cover
C -----------------------------------------------------------
C
          kgwanin(ji,jj) = kgwanin(ji,jj)*chemc(ji,jj,1)*1.e3* 
     &       (1-freeze(ji,jj))
         END DO
       END DO
C
C 2.5 COMPUTE GAS EXCHANGE COEFFICIENT FO O2 FROM LISS AND 
C      MERLIVAT EQUATIONS
C ---------------------------------------------------------
C
       DO jj = 1,jpj
         DO ji = 1,jpi
C
           ws = zwind(ji,jj)

           ttc = min(35.,tn(ji,jj,1))
           schmitt = 1953.4-128.0*ttc+3.9918*ttc**2
     &               -0.050091*ttc**3
C
           IF (ws.LE.3.6) THEN
               fugaci(ji,jj) = (a1*ws)*(schmitt/660.)**zexp1
           ENDIF
           IF ((ws.GT.3.6).AND.(ws.LE.13.)) THEN
               fugaci(ji,jj) = (a2*ws-b2)*(schmitt/660.)**zexp2
           ENDIF
           IF (ws.GT.13.) THEN
               fugaci(ji,jj) = (a3*ws-b3)*(schmitt/660.)**zexp2
           ENDIF
C
C CONVERT TO CM AND APPLY SEA ICE COVER
C -------------------------------------
C
           fugaci(ji,jj) = fugaci(ji,jj)/100./3600.*
     $         (1-freeze(ji,jj))*tmask(ji,jj,1)
C
#    if defined key_off_degrad
           fugaci(ji,jj) = exp(-rfact*fugaci(ji,jj)
     $         *facvol(ji,jj,1)/fse3t(ji,jj,1))
#    else
           fugaci(ji,jj) = exp(-rfact*fugaci(ji,jj)
     $         /fse3t(ji,jj,1))
#    endif
           
         ENDDO
       ENDDO
C

       DO jj = 1,jpj
         DO ji = 1,jpi
C
C Compute CO2 flux for the sea and air
C ------------------------------------
C
           alpco2 = chemc(ji,jj,1)
           fld = atcco2*tmask(ji,jj,1)*kgwanin(ji,jj)
           flu = h2co3(ji,jj)/alpco2
     &        *tmask(ji,jj,1)*kgwanin(ji,jj)

           tra(ji,jj,1,jpdic)= tra(ji,jj,1,jpdic)+(fld-flu)
     &        /1000./fse3t(ji,jj,1)
C
C Compute O2 flux 
C ---------------
C
          oxy16 = trn(ji,jj,1,jpoxy)
          flu16 = (-fugaci(ji,jj)+1)*fse3t(ji,jj,1)
     &            *(atcox*chemc(ji,jj,3)-oxy16)*
     &            tmask(ji,jj,1)/rfact 
          tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy)+flu16
     &            /fse3t(ji,jj,1)
 
C
C Save diagnostics
C ----------------
C
#    if defined key_trc_diaadd
          trc2d(ji,jj,1) = (fld-flu)
          trc2d(ji,jj,2) = flu16*1000.
          trc2d(ji,jj,3) = kgwanin(ji,jj)
          trc2d(ji,jj,4) = (fld-flu)/(kgwanin(ji,jj)+1.E-15)
#    endif
C
        END DO
      END DO
C
#endif
      RETURN
      END