source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trasbc.F90 @ 2117

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-09  (C. Ethe, G. Madec) Merge TRA-TRC
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   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_oce      ! ocean trends 
   USE trdtra          ! ocean trends
   USE in_out_manager  ! I/O manager
   USE prtctl          ! Print control
   USE sbcrnf          ! River runoff  
   USE sbcmod          ! ln_rnf  

   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. 
      !!
      !!      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 (Fwe=0), thus only Fwi remains:
      !!            sa = sa + emp * sn / e3t   for k=1
      !!         where emp, the surface freshwater budget (evaporation minus
      !!         precipitation minus runoff) given in kg/m2/s is divided
      !!         by 1035 kg/m3 (density of ocena 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 = 0, there is no surface flux of salt.
      !!
      !! ** 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')
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
      !!
      INTEGER  ::   ji, jj, jk      ! dummy loop indices  
      REAL(wp) ::   zta, zsa        ! local scalars, adjustment to temperature and salinity  
      REAL(wp) ::   zata, zasa      ! local scalars, calculations of automatic change to temp & sal due to vvl (done elsewhere)  
      REAL(wp) ::   zsrau, zse3t, zdep   ! local scalars, 1/density, 1/height of box, 1/height of effected water column  
      REAL(wp) ::   zdheat, zdsalt       ! total change of temperature and salinity  
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
      !!----------------------------------------------------------------------

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

      zsrau = 1. / rau0             ! initialization
#if defined key_zco
      zse3t = 1. / e3t_0(1)
#endif

      IF( l_trdtra )   THEN                    !* Save ta and sa trends
         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
      ENDIF

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

      ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff  
      DO jj = 2, jpj
         DO ji = fs_2, fs_jpim1   ! vector opt.
#if ! defined key_zco
            zse3t = 1. / fse3t(ji,jj,1)
#endif
            IF( lk_vvl) THEN
               ! temperature : heat flux and heat content of EMP flux
               zta = ( ro0cpr * qns(ji,jj) - emp(ji,jj) * zsrau * tsn(ji,jj,1,jp_tem) ) * zse3t
               ! Salinity    : concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration
               zsa = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t
            ELSE
               zta =  ro0cpr * qns(ji,jj) * zse3t                         ! temperature : heat flux 
               zsa =  emps(ji,jj) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t   ! salinity :  concent./dilut. effect 
            ENDIF
            tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta                  ! add the trend to the general tracer trend
            tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa
         END DO
      END DO

      IF( ln_rnf .AND. ln_rnf_att ) THEN        ! Concentration / dilution effect on (t,s) due to river runoff  
        DO jj = 1, jpj  
           DO ji = 1, jpi  
              rnf_dep(ji,jj) = 0.e0 
              DO jk = 1, rnf_mod_dep(ji,jj)                          ! recalculates rnf_dep to be the depth  
                rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk)    ! in metres to the bottom of the relevant grid box  
              END DO  
              zdep = 1. / rnf_dep(ji,jj)  
              zse3t= 1. / fse3t(ji,jj,1)  
              IF( rnf_tmp(ji,jj) == -999 )   rnf_tmp(ji,jj) = tsn(ji,jj,1,jp_tem)   ! if not specified set runoff temp to be sst  
  
              IF( rnf(ji,jj) > 0.e0 ) THEN  
  
                IF( lk_vvl ) THEN  
                  ! indirect flux, concentration or dilution effect : force a dilution effect in all levels 
                  zdheat = 0.e0  
                  zdsalt = 0.e0  
                  DO jk = 1, rnf_mod_dep(ji,jj)  
                    zta = -tsn(ji,jj,jk,jp_tem) * rnf(ji,jj) * zsrau * zdep  
                    zsa = -tsn(ji,jj,jk,jp_sal) * rnf(ji,jj) * zsrau * zdep  
                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend
                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
                    zdheat = zdheat + zta * fse3t(ji,jj,jk)  
                    zdsalt = zdsalt + zsa * fse3t(ji,jj,jk)  
                  END DO  
                  ! negate this total change in heat and salt content from top level    !!gm I don't understand this
                  zta = -zdheat * zse3t  
                  zsa = -zdsalt * zse3t  
                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta            ! add the trend to the general tracer trend
                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa
    
                  ! direct flux  
                  zta = rnf_tmp(ji,jj) * rnf(ji,jj) * zsrau * zdep  
                  zsa = rnf_sal(ji,jj) * rnf(ji,jj) * zsrau * zdep  
    
                  DO jk = 1, rnf_mod_dep(ji,jj)  
                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend
                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
                  END DO  
                ELSE  
                  DO jk = 1, rnf_mod_dep(ji,jj)  
                    zta = ( rnf_tmp(ji,jj) - tsn(ji,jj,jk,jp_tem) ) * rnf(ji,jj) * zsrau * zdep  
                    zsa = ( rnf_sal(ji,jj) - tsn(ji,jj,jk,jp_sal) ) * rnf(ji,jj) * zsrau * zdep  
                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend
                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
                  END DO  
                ENDIF  
  
              ELSEIF( rnf(ji,jj) < 0.e0) THEN   ! for use in baltic when flow is out of domain, want no change in temp and sal  
  
                IF( lk_vvl ) THEN  
                  ! calculate automatic adjustment to sal and temp due to dilution/concentraion effect   
                  zata = tsn(ji,jj,1,jp_tem) * rnf(ji,jj) * zsrau * zse3t  
                  zasa = tsn(ji,jj,1,jp_sal) * rnf(ji,jj) * zsrau * zse3t  
                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zata                  ! add the trend to the general tracer trend
                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zasa
                ENDIF  
  
              ENDIF  
  
           END DO  
        END DO  

      ELSE IF( ln_rnf ) THEN      ! Concentration dilution effect on (t,s) due to runoff without T, S and depth attributes


        DO jj = 2, jpj
           DO ji = fs_2, fs_jpim1   ! vector opt.
#if ! defined key_zco
              zse3t = 1. / fse3t(ji,jj,1)
#endif
              IF( lk_vvl) THEN
                 zta =    rnf(ji,jj) * zsrau * tsn(ji,jj,1,jp_tem) * zse3t       ! & cooling/heating effect of runoff
                 zsa =    0.e0                                                   ! No salinity concent./dilut. effect
              ELSE
                 zta =    0.0                                            ! temperature : heat flux 
                 zsa =  - rnf(ji,jj) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t     ! salinity :  concent./dilut. effect
              ENDIF
              tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta          ! add the trend to the general tracer trend
              tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa
           END DO
        END DO
 
      ENDIF  

      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_nsr, ztrdt )
         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_nsr, ztrds )
         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
      ENDIF
      !
      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
      !
   END SUBROUTINE tra_sbc

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