source: branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90 @ 1859

MODULE trasbc
   !!======================================================================
   !!                       ***  MODULE  trasbc  ***
   !! Ocean active tracers:  surface boundary condition
   !!======================================================================
   !! History :  OPA  !  1998-10  (G. Madec, G. Roullet, M. Imbard)  Original code
   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface
   !!  NEMO      1.0  !  2002-06  (G. Madec)  F90: Free form and module
   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   tra_sbc      : update the tracer trend at ocean surface
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and active tracers
   USE sbc_oce         ! surface boundary condition: ocean
   USE dom_oce         ! ocean space domain variables
   USE phycst          ! physical constant
   USE traqsr          ! solar radiation penetration
   USE trdmod          ! ocean trends 
   USE trdmod_oce      ! ocean variables trends
   USE in_out_manager  ! I/O manager
   USE prtctl          ! Print control

   IMPLICIT NONE
   PRIVATE

   PUBLIC   tra_sbc    ! routine called by step.F90

   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
   !! $Id$
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE tra_sbc ( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE tra_sbc  ***
      !!                   
      !! ** Purpose :   Compute the tracer surface boundary condition trend of
      !!      (flux through the interface, concentration/dilution effect)
      !!      and add it to the general trend of tracer equations.
      !!
      !! ** Method :
      !!      Following Roullet and Madec (2000), the air-sea flux can be divided 
      !!      into three effects: (1) Fext, external forcing; 
      !!      (2) Fwi, concentration/dilution effect due to water exchanged 
      !!         at the surface by evaporation, precipitations and runoff (E-P-R); 
      !!      (3) Fwe, tracer carried with the water that is exchanged. 
      !!            - salinity    : salt flux only due to freezing/melting
      !!            sa = sa +  fsalt / rau0 / e3t  for k=1
      !!
      !!      Fext, flux through the air-sea interface for temperature and salt: 
      !!            - temperature : heat flux q (w/m2). If penetrative solar
      !!         radiation q is only the non solar part of the heat flux, the
      !!         solar part is added in traqsr.F routine.
      !!            ta = ta + q /(rau0 rcp e3t)  for k=1
      !!            - salinity    : no salt flux
      !!
      !!      The formulation for Fwb and Fwi vary according to the free 
      !!      surface formulation (linear or variable volume). 
      !!      * Linear free surface
      !!            The surface freshwater flux modifies the ocean volume
      !!         and thus the concentration of a tracer and the temperature.
      !!         First order of the effect of surface freshwater exchange 
      !!         for salinity, it can be neglected on temperature (especially
      !!         as the temperature of precipitations and runoffs is usually
      !!         unknown).
      !!            - temperature : we assume that the temperature of both
      !!         precipitations and runoffs is equal to the SST, thus there
      !!         is no additional flux since in this case, the concentration
      !!         dilution effect is balanced by the net heat flux associated
      !!         to the freshwater exchange (Fwe+Fwi=0):
      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
      !!            - salinity    : evaporation, precipitation and runoff
      !!         water has a zero salinity  but there is a salt flux due to 
      !!         freezing/melting, thus:
      !!            sa = sa + emp * sn / rau0 / e3t   for k=1
      !!                    + fsalt    / rau0 / e3t
      !!         where emp, the surface freshwater budget (evaporation minus
      !!         precipitation minus runoff) given in kg/m2/s is divided
      !!         by rau0 = 1020 kg/m3 (density of sea water) to obtain m/s.    
      !!         Note: even though Fwe does not appear explicitly for 
      !!         temperature in this routine, the heat carried by the water
      !!         exchanged through the surface is part of the total heat flux
      !!         forcing and must be taken into account in the global heat
      !!         balance).
      !!      * nonlinear free surface (variable volume, lk_vvl)
      !!         contrary to the linear free surface case, Fwi is properly 
      !!         taken into account by using the true layer thicknesses to       
      !!         calculate tracer content and advection. There is no need to 
      !!         deal with it in this routine.
      !!           - temperature: Fwe=SST (P-E+R) is added to Fext.
      !!           - salinity:  Fwe sa = sa +  fsalt / rau0 / e3t.
      !!
      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
      !!                with the tracer surface boundary condition 
      !!              - save the trend it in ttrd ('key_trdtra')
      !!----------------------------------------------------------------------
      USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace   
      USE oce, ONLY :   ztrds => va   !  -  va     -      -   
      !!
      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
      !!
      INTEGER  ::   ji, jj        ! dummy loop indices
      REAL(wp) ::   z1_e3t_rau0   ! local scalars
      !!----------------------------------------------------------------------

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
         IF(lwp) WRITE(numout,*) '~~~~~~~ '
      ENDIF

      IF( l_trdtra ) THEN           ! Save ta and sa trends
         ztrdt(:,:,:) = ta(:,:,:) 
         ztrds(:,:,:) = sa(:,:,:) 
      ENDIF

!!gm useless staff ???
      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
!!gm

      
      IF( lk_vvl ) THEN          ! Variable Volume Layers case   ===>> heat content of mass flux in qns
         DO jj = 2, jpj
            DO ji = fs_2, fs_jpim1   ! vector opt.
               z1_e3t_rau0 = 1./ ( fse3t(ji,jj,1) * rau0 )
               ta(ji,jj,1) = ta(ji,jj,1) + z1_e3t_rau0 * qns (ji,jj)  * r1_rcp             ! non solar heat flux
               sa(ji,jj,1) = sa(ji,jj,1) + z1_e3t_rau0 * emps(ji,jj)                       ! salt flux (freezing/melting)
            END DO
         END DO
         !
      ELSE                       ! Constant Volume layers case   ===>> Concentration dillution effect
         DO jj = 2, jpj
            DO ji = fs_2, fs_jpim1   ! vector opt.
               z1_e3t_rau0 = 1./ ( fse3t(ji,jj,1) * rau0 )
               ta(ji,jj,1) = ta(ji,jj,1) + z1_e3t_rau0 * (  qns (ji,jj) * r1_rcp       &   ! non solar heat flux
                  &                                       + emp (ji,jj) * tn(ji,jj,1)  )   ! concent./dilut. effect
               sa(ji,jj,1) = sa(ji,jj,1) + z1_e3t_rau0 * (  emps(ji,jj)                &   ! salt flux (freezing/melting)
                  &                                       + emp (ji,jj) * sn(ji,jj,1)  )   ! concent./dilut. effect
            END DO
         END DO
      ENDIF

      IF( l_trdtra ) THEN      ! save the sbc trends for diagnostic
         ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
         ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)
         CALL trd_mod(ztrdt, ztrds, jptra_trd_nsr, 'TRA', kt)
      ENDIF
      !
      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' sbc  - Ta: ', mask1=tmask,   &
         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
      !
   END SUBROUTINE tra_sbc

   !!======================================================================
END MODULE trasbc