New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 1859 – NEMO

Changeset 1859


Ignore:
Timestamp:
2010-05-06T10:40:07+02:00 (14 years ago)
Author:
gm
Message:

ticket:#665 step 2 & 3: heat content in qns & new forcing terms

Location:
branches/DEV_r1837_mass_heat_salt_fluxes/NEMO
Files:
17 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90

    r1858 r1859  
    1515   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    1616   !!---------------------------------------------------------------------- 
    17    !! * Modules used 
    1817   USE par_ice_2          ! LIM sea-ice parameters 
    1918 
     
    3433   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation 
    3534   REAL(wp), PUBLIC ::   rdt_ice            !: ice time step 
     35   REAL(wp), PUBLIC ::   r1_rdt_ice         !: =1/rdt_ice 
    3636   REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic 
    3737   REAL(wp), PUBLIC ::   alpha  = 0.5       !: coefficient for semi-implicit coriolis 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/iceini_2.F90

    r1857 r1859  
    1111   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
     13   !!   ice_init_2     : LIM-2 sea-ice model initialization 
    1314   !!---------------------------------------------------------------------- 
    14    !!   ice_init_2       : sea-ice model initialization 
    15    !!   ice_run_2        : Definition some run parameter for ice model 
    16    !!---------------------------------------------------------------------- 
    17    USE dom_oce 
    18    USE dom_ice_2 
    19    USE sbc_oce         ! surface boundary condition: ocean 
    20    USE sbc_ice         ! surface boundary condition: ice 
    21    USE phycst          ! Define parameters for the routines 
    22    USE ice_2 
    23    USE limmsh_2 
    24    USE limistate_2 
    25    USE limrst_2    
    26    USE in_out_manager 
     15   USE dom_oce         ! ocean domain 
     16   USE dom_ice_2       ! LIM-2 : ice domain 
     17   USE sbc_oce         ! ocean surface boundary condition 
     18   USE sbc_ice         ! ice   surface boundary condition 
     19   USE phycst          ! physical constant 
     20   USE ice_2           ! LIM-2 variables 
     21   USE limmsh_2        ! LIM-2 mesh 
     22   USE limistate_2     ! LIM-2 inital state 
     23   USE limrst_2        ! LIM-2 restart 
     24   USE in_out_manager  ! I/O manager 
    2725       
    2826   IMPLICIT NONE 
    2927   PRIVATE 
    3028 
    31    PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90 
     29   PUBLIC   ice_init_2   ! called by sbcice_lim_2.F90 
    3230 
    3331   !!---------------------------------------------------------------------- 
    34    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
     32   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
    3533   !! $Id$  
    36    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     34   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3735   !!---------------------------------------------------------------------- 
    3836 
     
    4341      !!                  ***  ROUTINE ice_init_2  *** 
    4442      !! 
    45       !! ** purpose :    
     43      !! ** purpose :   LIM-2 sea-ice model initialisation 
     44      !! 
     45      !! ** input   :   namelist_ice file 
     46      !!                namelist namicerun (inside namelist_ice file) 
    4647      !!---------------------------------------------------------------------- 
    47       ! 
    48       ! Open the namelist file  
    49       CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )       
    50       CALL ice_run_2                    !  read in namelist some run parameters 
    51                   
    52       ! Louvain la Neuve Ice model 
    53       rdt_ice = nn_fsbc * rdttra(1) 
    54  
    55       CALL lim_msh_2                  ! ice mesh initialization 
    56       
    57       ! Initial sea-ice state 
    58       IF( .NOT.ln_rstart ) THEN 
    59          CALL lim_istate_2            ! start from rest: sea-ice deduced from sst 
    60       ELSE 
    61          CALL lim_rst_read_2          ! start from a restart file 
    62       ENDIF 
    63        
    64       tn_ice(:,:,1) = sist(:,:)         ! initialisation of ice temperature    
    65       fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction     
    66       ! 
    67    END SUBROUTINE ice_init_2 
    68  
    69  
    70    SUBROUTINE ice_run_2 
    71       !!------------------------------------------------------------------- 
    72       !!                  ***  ROUTINE ice_run_2 *** 
    73       !!                  
    74       !! ** Purpose :   Definition some run parameter for ice model 
    75       !! 
    76       !! ** Method  :   Read the namicerun namelist and check the parameter  
    77       !!       values called at the first timestep (nit000) 
    78       !! 
    79       !! ** input   :   Namelist namicerun 
    80       !!------------------------------------------------------------------- 
    8148      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, ln_limdmp, acrit, hsndif, hicdif 
    8249      !!------------------------------------------------------------------- 
    83       !                     
    84       REWIND ( numnam_ice )                       ! Read Namelist namicerun  
     50      ! 
     51      !                                    ! Open the ice namelist file  
     52      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )    
     53      ! 
     54      REWIND ( numnam_ice )                ! Read Namelist namicerun  
    8555      READ   ( numnam_ice , namicerun ) 
    86  
    87       IF(lwp) THEN 
     56      ! 
     57      IF(lwp) THEN                         ! control print 
    8858         WRITE(numout,*) 
    8959         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
     
    9565         WRITE(numout,*) '   computation of temp. in ice  (=0) or not (=9999) hicdif = ', hicdif 
    9666      ENDIF 
     67      !          
     68      rdt_ice    = nn_fsbc * rdttra(1)     ! set ice time step to surface tracer time step 
     69      r1_rdt_ice = 1.e0 / rdt_ice 
    9770      ! 
    98    END SUBROUTINE ice_run_2 
     71      CALL lim_msh_2                       ! ice mesh initialization 
     72      ! 
     73      !                                    ! Initial sea-ice state 
     74      IF( .NOT.ln_rstart ) THEN   ;   CALL lim_istate_2        ! start from rest: sea-ice deduced from sst 
     75      ELSE                        ;   CALL lim_rst_read_2      ! start from a restart file 
     76      ENDIF 
     77      ! 
     78      tn_ice(:,:,1) = sist(:,:)            ! set initial ice temperature   
     79      ! 
     80      fr_i  (:,:) = 1.0 - frld(:,:)        ! set initial ice fraction     
     81      ! 
     82   END SUBROUTINE ice_init_2 
    9983 
    10084#else 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90

    r1858 r1859  
    44   !!           computation of the flux at the sea ice/ocean interface 
    55   !!====================================================================== 
    6    !! History : LIM  ! 2000-01 (H. Goosse) Original code 
    7    !!           2.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90 
    8    !!            -   ! 2006-07 (G. Madec) surface module 
    9    !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
     6   !! History :  1.0  !  2000-01  (H. Goosse) Original code 
     7   !!            2.0  !  2002-07  (C. Ethe, G. Madec) re-writing F90 
     8   !!             -   !  2006-07  (G. Madec) surface module 
     9   !!             -   !  2008-07  (C. Talandier,G.  Madec) 2D fields for soce and sice 
     10   !!            2.1  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_lim2 
     
    1718   USE par_oce          ! ocean parameters 
    1819   USE dom_oce          ! ocean domain 
    19    USE sbc_ice          ! surface boundary condition 
    20    USE sbc_oce          ! surface boundary condition 
     20   USE sbc_ice          ! ice   surface boundary condition 
     21   USE sbc_oce          ! ocean surface boundary condition 
    2122   USE phycst           ! physical constants 
    22    USE ice_2            ! LIM sea-ice variables 
     23   USE albedo           ! albedo parameters 
     24   USE ice_2            ! LIM-2 sea-ice variables 
    2325 
    2426   USE lbclnk           ! ocean lateral boundary condition 
    2527   USE in_out_manager   ! I/O manager 
     28   USE iom              !  
     29   USE prtctl           ! Print control 
    2630   USE diaar5, ONLY :   lk_diaar5 
    27    USE iom              !  
    28    USE albedo           ! albedo parameters 
    29    USE prtctl           ! Print control 
    3031   USE cpl_oasis3, ONLY : lk_cpl 
    3132 
     
    3334   PRIVATE 
    3435 
    35    PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2 
     36   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2 
    3637 
    3738   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values 
    3839   REAL(wp)  ::   rzero  = 0.e0     
    3940   REAL(wp)  ::   rone   = 1.e0 
    40    REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r 
    41    REAL(wp), DIMENSION(jpi,jpj)  ::   sice_r 
     41   REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r, sice_r   ! ocean and ice 2D constant salinity fields (used if lk_vvl=F) 
    4242 
    4343   !! * Substitutions 
    4444#  include "vectopt_loop_substitute.h90" 
    4545   !!---------------------------------------------------------------------- 
    46    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     46   !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
    4747   !! $Id$ 
    4848   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    6262      !!              - Update  
    6363      !!      
    64       !! ** Outputs : - qsr     : sea heat flux:     solar  
    65       !!              - qns     : sea heat flux: non solar 
    66       !!              - emp     : freshwater budget: volume flux  
    67       !!              - emps    : freshwater budget: concentration/dillution  
     64      !! ** Outputs : - qsr     : solar heat flux 
     65      !!              - qns     : non-solar heat flux including heat content of mass flux 
     66      !!              - emp     : mass flux 
     67      !!              - emps    : salt flux due to Freezing/Melting  
    6868      !!              - utau    : sea surface i-stress (ocean referential) 
    6969      !!              - vtau    : sea surface j-stress (ocean referential) 
     
    7575      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 
    7676      !!--------------------------------------------------------------------- 
    77       INTEGER ::   kt    ! number of iteration 
     77      INTEGER, INTENT(in) ::   kt    ! number of iteration 
    7878      !! 
    79       INTEGER  ::   ji, jj           ! dummy loop indices 
    80       INTEGER  ::   ifvt, i1mfr, idfr               ! some switches 
    81       INTEGER  ::   iflt, ial, iadv, ifral, ifrdv 
    82       INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers 
    83       REAL(wp) ::   zrdtir           ! 1. / rdt_ice 
    84       REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux 
    85       REAL(wp) ::   zinda            ! switch for testing the values of ice concentration 
    86       REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface 
    87       REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface 
    88       REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points 
    89       REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points 
    90       REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity 
    91 ! interface 2D --> 3D 
    92       REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb     ! albedo of ice under overcast sky 
    93       REAL(wp), DIMENSION(jpi,jpj,1) ::   zalbp    ! albedo of ice under clear sky 
     79      INTEGER  ::   ji, jj                     ! dummy loop indices 
     80      INTEGER  ::   ifvt, idfr , iadv, i1mfr   ! local integers 
     81      INTEGER  ::   iflt, ifrdv, ial , ifral   !   -      - 
     82      INTEGER  ::   ii0, ii1, ij0, ij1         !   -      - 
     83      REAL(wp) ::   zqsr, zqns, zqhc, zemp     ! local scalars 
     84      REAL(wp) ::   zinda, zswitch, zcd        !   -      - 
     85      REAL(wp) ::   zfrldu, zutau, zu_io       !   -      - 
     86      REAL(wp) ::   zfrldv, zvtau, zv_io       !   -      - 
     87      REAL(wp) ::   zemp_snw, zfmm, zfsalt     !   -      - 
    9488      REAL(wp) ::   zsang, zmod, zztmp, zfm 
    95       REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! component of ocean stress below sea-ice at I-point 
    96       REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi           ! module    of ocean stress below sea-ice at I-point 
    97       REAL(wp), DIMENSION(jpi,jpj) ::   zqnsoce          ! save qns before its modification by ice model 
    98  
     89      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp    ! 3D workspace 
     90      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! 2D workspace 
     91      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi, zqnsoce  !  -      - 
    9992      !!--------------------------------------------------------------------- 
    100       
    101       zrdtir = 1. / rdt_ice 
    102        
     93            
    10394      IF( kt == nit000 ) THEN 
    10495         IF(lwp) WRITE(numout,*) 
    10596         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition' 
    10697         IF(lwp) WRITE(numout,*) '~~~~~~~~~   ' 
    107  
    108          soce_r(:,:) = soce 
    109          sice_r(:,:) = sice 
    110          ! 
    111          IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN 
    112             !                                        ! ======================= 
    113             !                                        !  ORCA_R2 configuration 
    114             !                                        ! ======================= 
    115             ii0 = 145   ;   ii1 = 180        ! Baltic Sea 
     98         !                              ! 2D fields for constant ice and ocean salinities 
     99         soce_r(:,:) = soce             !    in order to use different value in the Baltic sea 
     100         sice_r(:,:) = sice             !    which is much less salty than polar regions 
     101         ! 
     102         IF( cp_cfg == "orca" ) THEN    ! ORCA configuration 
     103            IF( jp_cfg == 2       ) THEN     !  ORCA_R2 configuration 
     104            ii0 = 145   ;   ii1 = 180              ! Baltic Sea 
    116105            ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 
    117                                             sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 2.e0 
     106                                            sice_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 
     107!!gm add here the R1 R05 and R025  cases 
     108!!          ELSEIF( jp_cfg == 1   ) THEN           !  ORCA_R1   configuration 
     109!!          ELSEIF( jp_cfg == 05  ) THEN           !  ORCA_R05  configuration 
     110!!          ELSEIF( jp_cfg == 025 ) THEN           !  ORCA_R025 configuration 
     111!! 
     112!!gm or better introduce the baltic change as a function of lat/lon of the baltic sea 
     113!!end gm 
     114            ENDIF 
    118115         ENDIF 
    119116         ! 
    120117      ENDIF 
    121118 
    122       !------------------------------------------! 
    123       !      heat flux at the ocean surface      ! 
    124       !------------------------------------------! 
    125  
    126 !!gm 
    127 !!gm CAUTION    
    128 !!gm re-verifies the non solar expression, especially over open ocen 
    129 !!gm 
    130       zqnsoce(:,:) = qns(:,:) 
     119      zqnsoce(:,:) = qns(:,:)      ! save non-solar flux prior to its modification by ice-ocean fluxes (diag.) 
     120      ! 
     121      zswitch = 1       ! standard levitating sea-ice : salt exchange only 
     122      ! 
     123!!gm ice embedment 
     124!      SELECT CASE( nn_ice_embd )       ! levitating/embedded sea-ice 
     125!      CASE( 0    )   ;   zswitch = 1       ! standard levitating sea-ice : salt exchange only 
     126!      CASE( 1, 2 )   ;   zswitch = 0       ! other levitating sea-ice or embedded sea-ice : salt and volume fluxes 
     127!      END SELECT                           !     
     128!!gm end embedment 
     129 
    131130      DO jj = 1, jpj 
    132131         DO ji = 1, jpi 
     132            !                          !------------------------------------------! 
     133            !                          !      heat flux at the ocean surface      ! 
     134            !                          !------------------------------------------! 
    133135            zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) ) 
    134136            ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) ) 
     
    138140            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr 
    139141            iadv    = ( 1  - i1mfr ) * zinda 
    140             ifral   = ( 1  - i1mfr * ( 1 - ial ) )    
    141             ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv  
    142  
    143 !!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time) 
    144 !!$ 
    145 !!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time) 
    146 !!$ 
    147 !!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ??? 
    148 !!$            ELSE                             ;   ifvt = 0. 
    149 !!$            ENDIF 
    150 !!$ 
    151 !!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current 
    152 !!$            ELSE                                     ;   idfr = 1.    
    153 !!$            ENDIF 
    154 !!$ 
    155 !!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current 
    156 !!$ 
     142            ifral   = ( 1  - i1mfr * ( 1 - ial ) ) 
     143            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
     144 
     145!!gm  attempt to understand and comment the tricky flags used here.... 
     146! 
     147!gm      zinda   = 1.0 - AINT( pfrld(ji,jj) )    ! = 0. free-ice ocean else 1. (after ice adv, but before ice thermo) 
     148!gm      i1mfr   = 1.0 - AINT(  frld(ji,jj) )    ! = 0. free-ice ocean else 1. (after ice thermo) 
     149! 
     150!gm      IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda   ! = 1. if (snow and no ice at previous time) else 0. ??? 
     151!gm      ELSE                             ;   ifvt = 0.      ! correspond to a overmelting of snow in surface ablation 
     152!gm      ENDIF                                               !  
     153! 
     154!gm      IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermo 
     155!gm      ELSE                                     ;   idfr = 1.    
     156!gm      ENDIF 
     157! 
     158!!$      iflt    = zinda  * (1 - i1mfr) * (1 - ifvt ) ! = 1. if ice (not only snow) at previous and pure ocean at current 
     159! 
    157160!!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr 
    158161!!$!                 snow no ice   ice         ice or nothing  lead fraction increases 
    159162!!$!                 at previous   now           at previous 
    160163!!$!                -> ice aera increases  ???         -> ice aera decreases ??? 
    161 !!$ 
     164! 
    162165!!$            iadv    = ( 1  - i1mfr ) * zinda   
    163166!!$!                     pure ocean      ice at 
    164167!!$!                     at current      previous 
    165168!!$!                        -> = 1. if ice disapear between previous and current 
    166 !!$ 
     169! 
    167170!!$            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
    168171!!$!                            ice at     ??? 
    169172!!$!                            current          
    170173!!$!                         -> ??? 
    171 !!$  
     174! 
    172175!!$            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv  
    173176!!$!                                                    ice disapear                            
    174 !!$ 
    175 !!$ 
    176  
    177             !   computation the solar flux at ocean surface 
     177! 
     178            ! 
     179            ! - computation the solar flux at ocean surface 
    178180#if defined key_coupled  
    179181            zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) ) 
     
    181183            zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj) 
    182184#endif             
    183             !  computation the non solar heat flux at ocean surface 
     185            ! 
     186            ! - computation the non solar heat flux at ocean surface 
    184187            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads 
    185                &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            & 
    186                &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    & 
    187                &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir 
    188  
    189             fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ??? 
    190              
    191             qsr  (ji,jj) = zqsr                                          ! solar heat flux  
    192             qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux 
     188               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                                & 
     189               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdt_ice    & 
     190               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdt_ice 
     191 
     192            ! - store residual heat flux (put in the ocean at the next time-step) 
     193            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)   ! ??? 
     194            ! 
     195            ! - heat content of mass exchanged between ocean and sea-ice 
     196            zqhc = ( rdq_snw(ji,jj) + rdq_ice(ji,jj) ) * r1_rdt_ice    ! heat flux due to sown & ice heat content exchanges 
     197            !             
     198            qsr(ji,jj) = zqsr                                          ! solar heat flux  
     199            qns(ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux 
     200   
     201            !                          !------------------------------------------! 
     202            !                          !      mass flux at the ocean surface      ! 
     203            !                          !------------------------------------------! 
     204            ! 
     205            ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area) 
     206#if defined key_coupled 
     207            !                                                       ! coupled mode:  
     208            zemp = + emp_tot(ji,jj)                              &       ! net mass flux over the grid cell (ice+ocean area) 
     209               &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )              ! minus the mass flux intercepted by sea-ice 
     210#else 
     211            !                                                       ! forced  mode:  
     212            zemp = + emp(ji,jj)     *         frld(ji,jj)      &         ! mass flux over open ocean fraction  
     213               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &         ! liquid precip. over ice reaches directly the ocean 
     214               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &         ! snow is intercepted by sea-ice (previous frld) 
     215#endif             
     216            ! 
     217            ! mass flux at the ocean/ice interface (sea ice fraction) 
     218            zemp_snw = rdm_snw(ji,jj) * r1_rdt_ice                  ! snow melting = pure water that enters the ocean 
     219            zfmm     = rdm_ice(ji,jj) * r1_rdt_ice                  ! Freezing minus Melting (F-M) 
     220 
     221            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s] 
     222            zfsalt = - sice_r(ji,jj) * zfmm                         ! F-M salt exchange 
     223            zcd    =   soce_r(ji,jj) * zfmm                         ! concentration/dilution term due to F-M 
     224            ! 
     225            ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux 
     226            ! salt and mass fluxes : non concentartion dilution term in salt flux  and add F-M term in volume flux 
     227            emps(ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange) 
     228            emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (- F/M mass flux if no ice/ocean mass exchange) 
     229            ! 
    193230         END DO 
    194231      END DO 
     232 
    195233 
    196234      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )       
    197235      CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )       
    198236      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) ) 
    199  
    200       !------------------------------------------! 
    201       !      mass flux at the ocean surface      ! 
    202       !------------------------------------------! 
    203  
    204 !!gm 
    205 !!gm CAUTION    
    206 !!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm 
    207 !!gm 
    208       DO jj = 1, jpj 
    209          DO ji = 1, jpi 
    210              
    211 #if defined key_coupled 
    212           zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  
    213              &   + rdm_snw(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting  
    214 #else 
    215 !!$            !  computing freshwater exchanges at the ice/ocean interface 
    216 !!$            zpme = - evap(ji,jj)    *   frld(ji,jj)           &   !  evaporation over oceanic fraction 
    217 !!$               &   + tprecip(ji,jj)                           &   !  total precipitation 
    218 !!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice 
    219 !!$               &   - rdm_snw(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting  
    220             !  computing freshwater exchanges at the ice/ocean interface 
    221             zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction  
    222                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean 
    223                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step 
    224                &   + rdm_snw(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting  
    225                !                                                   !  ice-covered fraction: 
    226 #endif             
    227  
    228             !  computing salt exchanges at the ice/ocean interface 
    229             zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdm_ice(ji,jj) * zrdtir )  
    230              
    231             !  converting the salt flux from ice to a freshwater flux from ocean 
    232             zfm  = zfons / ( sss_m(ji,jj) + epsi16 ) 
    233              
    234             emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution) 
    235             emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution) 
    236  
    237          END DO 
    238       END DO 
    239237 
    240238      IF( lk_diaar5 ) THEN 
     
    250248      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case) 
    251249         !                                         ! otherwise the atmosphere-ocean stress is used everywhere 
    252  
     250         ! 
    253251         ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point) 
    254252!CDIR NOVERRCHK 
     
    290288            END DO 
    291289         END DO 
    292  
    293          ! boundary condition on the stress (utau,vtau,taum) 
    294          CALL lbc_lnk( utau, 'U', -1. ) 
    295          CALL lbc_lnk( vtau, 'V', -1. ) 
     290         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )         ! lateral boundary conditions 
    296291         CALL lbc_lnk( taum, 'T',  1. ) 
    297  
     292         ! 
    298293      ENDIF 
    299294 
    300295      !-----------------------------------------------! 
    301       !   Coupling variables                          ! 
     296      !   Storing the transmitted variables           ! 
    302297      !-----------------------------------------------! 
    303  
    304       IF ( lk_cpl ) THEN            
    305          ! Ice surface temperature  
     298!!gm where this is done ?????   ==>>> limthd_2   not logic ??? 
     299!!gm      fr_i(:,:) = 1.0 - frld(:,:)       ! sea-ice fraction 
     300!!gm 
     301 
     302      IF ( lk_cpl ) THEN      ! coupled mode : 
    306303         tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature        
    307          ! Computation of snow/ice and ocean albedo 
     304         !                                  ! snow/ice and ocean albedo 
    308305         CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb ) 
    309306         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys) 
     307         ! 
    310308         CALL iom_put( "icealb_cea", alb_ice(:,:,1) * fr_i(:,:) )  ! ice albedo 
    311309      ENDIF 
     
    318316         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ') 
    319317      ENDIF  
    320     
    321     END SUBROUTINE lim_sbc_2 
     318      ! 
     319   END SUBROUTINE lim_sbc_2 
    322320 
    323321#else 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90

    r1858 r1859  
    88   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
    99   !!---------------------------------------------------------------------- 
    10    !!   LIM 2.0, UCL-LOCEAN-IPSL (2005) 
     10   !! NEMO/LIM 3.3, UCL-LOCEAN-IPSL (2010) 
    1111   !! $Id$ 
    12    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
     12   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    1313   !!---------------------------------------------------------------------- 
    14    !! * Modules used 
    1514   USE par_ice_2 
    1615 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/ice.F90

    r1471 r1859  
    11MODULE ice 
     2   !!====================================================================== 
     3   !!                        ***  MODULE ice  *** 
     4   !! LIM-3 sea ice :   variables defined in memory 
     5   !!===================================================================== 
     6   !! History :  3.0  ! 2008-03  (M. Vancoppenolle)  LIM3 original code 
     7   !!            3.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
     8   !!---------------------------------------------------------------------- 
    29#if defined key_lim3 
    310   !!---------------------------------------------------------------------- 
    411   !!   'key_lim3' :                                   LIM3 sea-ice model 
    512   !!---------------------------------------------------------------------- 
    6    !! History : 
    7    !!   2.0  !  03-08  (C. Ethe)  F90: Free form and module 
    8    !!   3.0  !  08-03  (M. Vancoppenolle) : LIM3 ! 
    9    !!---------------------------------------------------------------------- 
    10    !!  LIM 3.0, UCL-LOCEAN-IPSL (2005) 
    11    !! $Id$ 
    12    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    13    !!---------------------------------------------------------------------- 
    14    !! * Modules used 
    1513   USE par_ice          ! LIM sea-ice parameters 
    1614 
    1715   IMPLICIT NONE 
    1816   PRIVATE 
    19    !! 
    20    !!====================================================================== 
    21    !!                        ***  MODULE ice  *** 
    22    !! 
    23    !!                           ************** 
    24    !!                           * L I M  3.0 * 
    25    !!                           ************** 
    26    !! 
    27    !!                         ''in ice we trust'' 
    28    !! 
    29    !!                   This module contains the sea ice  
    30    !!                 diagnostics variables of ice defined  
    31    !!                             in memory 
    32    !! 
    33    !!====================================================================== 
    34    !! 
     17 
     18   !!---------------------------------------------------------------------- 
    3519   !! LIM3 by the use of sweat, agile fingers and sometimes brain juice,  
    3620   !!  was developed in Louvain-la-Neuve by :  
    37    !! 
    3821   !!    * Martin Vancoppenolle (UCL-ASTR, Belgium) 
    3922   !!    * Sylvain Bouillon (UCL-ASTR, Belgium) 
     
    4124   !!  
    4225   !! Based on extremely valuable earlier work by 
    43    !! 
    4426   !!    * Thierry Fichefet 
    4527   !!    * Hugues Goosse 
    4628   !! 
    4729   !! The following persons also contributed to the code in various ways 
    48    !! 
    4930   !!    * Gurvan Madec, Claude Talandier, Christian Ethe  
    5031   !!      and Rachid Benshila (LOCEAN-IPSL, France) 
    5132   !!    * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany) 
    52    !!    * Bill Lipscomb (LANL), Cecilia Bitz (UWa)  
    53    !!      and Elisabeth Hunke (LANL), USA. 
     33   !!    * Bill Lipscomb (LANL), Cecilia Bitz (UWa) and Elisabeth Hunke (LANL), USA. 
    5434   !!  
    5535   !! (c) UCL-ASTR, 2005-2008 
    5636   !! 
    57    !! For more info, the interested user is kindly invited to consult the 
    58    !! following references 
     37   !! For more info, the interested user is kindly invited to consult the following references 
    5938   !!    For model description and validation : 
    6039   !!    * Vancoppenolle et al., Ocean Modelling, 2008a. 
    6140   !!    * Vancoppenolle et al., Ocean Modelling, 2008b. 
    62    !! 
    6341   !!    For a specific description of EVP : 
    6442   !!    * Bouillon et al., in prep for 2008. 
    65    !! 
    66    !!    Or the reference manual, that should be available by 2009 
    6743   !! 
    6844   !!====================================================================== 
     
    183159   !!===================================================================== 
    184160 
    185    LOGICAL, PUBLIC ::    & 
    186       con_i = .false.           ! switch for conservation test 
     161   LOGICAL, PUBLIC ::   con_i = .false.   ! switch for conservation test 
    187162 
    188163   !!-------------------------------------------------------------------------- 
     
    310285      phicif ,   &  !: Old ice thickness 
    311286      fbif   ,   &  !: Heat flux at the ice base 
    312       rdmsnif,   &  !: Variation of snow mass 
    313       rdmicif,   &  !: Variation of ice mass 
     287      rdm_snw,   &  !: Variation of snow mass over 1 time step           [Kg/m2] 
     288      rdq_snw,   &  !: heat content associated to rdm_snw                [J/m2] 
     289      rdm_ice,   &  !: Variation of ice  mass over 1 time step           [Kg/m2] 
     290      rdq_ice,   &  !: heat content associated to rdm_ice                [J/m2] 
    314291      qldif  ,   &  !: heat balance of the lead (or of the open ocean) 
    315292      qcmif  ,   &  !: Energy needed to bring the ocean surface layer until its freezing  
     
    506483#endif 
    507484 
     485   !!---------------------------------------------------------------------- 
     486   !! NEMO/LIM 3.0, UCL-LOCEAN-IPSL (2010) 
     487   !! $Id$ 
     488   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    508489   !!====================================================================== 
    509490END MODULE ice 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/thd_ice.F90

    r1465 r1859  
    22   !!====================================================================== 
    33   !!                       ***  MODULE thd_ice  *** 
    4    !! LIM sea-ice :   Ice thermodynamics in 1D 
     4   !! LIM-3 sea-ice :   Ice thermodynamics in 1D 
    55   !!===================================================================== 
    6    !! History : 
    7    !!   2.0  !  02-11  (C. Ethe)  F90: Free form and module 
     6   !! History :  3.0  ! 2008-03  (M. Vancoppenolle)  LIM3 original code 
     7   !!            3.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
    88   !!---------------------------------------------------------------------- 
    9    !!   LIM 3.0, UCL-LOCEAN-IPSL (2008) 
    10    !! $Id$ 
    11    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    12    !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    149   USE par_ice 
    1510 
     
    7368      at_i_b      ,     &  !:    "                  "      frld 
    7469      fbif_1d     ,     &  !:    "                  "      fbif 
    75       rdmicif_1d  ,     &  !:    "                  "      rdmicif 
    76       rdmsnif_1d  ,     &  !:    "                  "      rdmsnif 
     70      rdm_ice_1d  ,     &  !:    "                  "      rdm_ice 
     71      rdq_ice_1d  ,     &  !:    "                  "      rdm_ice 
     72      rdm_snw_1d  ,     &  !:    "                  "      rdm_snw 
     73      rdq_snw_1d  ,     &  !:    "                  "      rdm_snw 
    7774      qlbbq_1d    ,     &  !:    "                  "      qlbsbq 
    7875      dmgwi_1d    ,     &  !:    "                  "      dmgwi 
     
    149146      ftotal_fin               !: final total heat flux 
    150147 
    151    REAL(wp), PUBLIC, DIMENSION(jpij,0:nlay_s) ::   &  !: 
    152       fc_s 
    153    REAL(wp), PUBLIC, DIMENSION(jpij,0:jkmax)  ::   &  !: 
    154       fc_i 
    155    REAL(wp), PUBLIC, DIMENSION(jpij,nlay_s) ::   &  !: 
    156       de_s_lay 
    157    REAL(wp), PUBLIC, DIMENSION(jpij,jkmax)  ::   &  !: 
    158       de_i_lay 
    159    INTEGER , PUBLIC ::                           & 
    160       jiindex_1d   ! 1D index of debugging point 
     148   REAL(wp), PUBLIC, DIMENSION(jpij,0:nlay_s) ::   fc_s 
     149   REAL(wp), PUBLIC, DIMENSION(jpij,0:jkmax)  ::   fc_i 
     150   REAL(wp), PUBLIC, DIMENSION(jpij,nlay_s)   ::   de_s_lay 
     151   REAL(wp), PUBLIC, DIMENSION(jpij,jkmax)    ::   de_i_lay 
     152   INTEGER , PUBLIC                           ::   jiindex_1d   ! 1D index of debugging point 
    161153 
     154   !!---------------------------------------------------------------------- 
     155   !! NEMO/LIM 3.0, UCL-LOCEAN-IPSL (2010) 
     156   !! $Id$ 
     157   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    162158   !!====================================================================== 
    163159END MODULE thd_ice 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcana.F90

    r1732 r1859  
    66   !! History :  3.0   ! 2006-06  (G. Madec)  Original code 
    77   !!            3.2   ! 2009-07  (G. Madec)  Style only 
     8   !!            3.3  !  2010-07  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    89   !!---------------------------------------------------------------------- 
    910 
     
    3940#  include "vectopt_loop_substitute.h90" 
    4041   !!---------------------------------------------------------------------- 
    41    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     42   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    4243   !! $Id$ 
    4344   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    6061      !! 
    6162      !! ** Action  : - set the ocean surface boundary condition, i.e.   
    62       !!                   utau, vtau, taum, wndm, qns, qsr, emp, emps 
     63      !!                   utau, vtau, taum, wndm, qns, qsr, emp 
    6364      !!---------------------------------------------------------------------- 
    6465      INTEGER, INTENT(in) ::   kt       ! ocean time step 
     
    8889         ! 
    8990         nn_tau000 = MAX( nn_tau000, 1 )   ! must be >= 1 
    90          qns   (:,:) = rn_qns0 
    91          qsr   (:,:) = rn_qsr0 
    92          emp   (:,:) = rn_emp0 
    93          emps  (:,:) = rn_emp0 
     91         emp(:,:) = rn_emp0 
     92         qns(:,:) = rn_qns0 - emp(:,:) * sst_m(:,:) * rcp      ! including heat content associated with mass flux at SST 
     93         qsr(:,:) = rn_qsr0 
    9494         ! 
    9595      ENDIF 
     
    123123      !! 
    124124      !! ** Action  : - set the ocean surface boundary condition, i.e.    
    125       !!                   utau, vtau, taum, wndm, qns, qsr, emp, emps 
     125      !!                   utau, vtau, taum, wndm, qns, qsr, emp 
    126126      !! 
    127127      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000. 
     
    204204         END DO 
    205205      END DO 
    206       emps(:,:) = emp(:,:) 
    207206 
    208207      ! Compute the emp flux such as its integration on the whole domain at each time is zero 
     
    226225      ENDIF 
    227226 
    228       !salinity terms 
    229       emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1) 
    230       emps(:,:) = emp(:,:) 
     227       
     228      ! freshwater (mass flux) and update of qns with heat content of emp 
     229      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average) 
     230      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST 
    231231 
    232232 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

    r1732 r1859  
    99   !!            3.0  !  2008-03 (C. Talandier, G. Madec) surface module + LIM3 
    1010   !!            3.2  !  2009-04 (B. Lemaire) Introduce iom_put 
     11   !!            3.3  !  2010-05 (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    1112   !!---------------------------------------------------------------------- 
    1213 
    1314   !!---------------------------------------------------------------------- 
    14    !!   sbc_blk_clio   : CLIO bulk formulation: read and update required input fields 
    15    !!   blk_clio_oce   : ocean CLIO bulk formulea: compute momentum, heat and freswater fluxes for the ocean 
    16    !!   blk_ice_clio   : ice   CLIO bulk formulea: compute momentum, heat and freswater fluxes for the sea-ice 
     15   !!   sbc_blk_clio     : CLIO bulk formulation: read and update required input fields 
     16   !!   blk_clio_oce     : ocean CLIO bulk formulea: compute momentum, heat and freswater fluxes for the ocean 
     17   !!   blk_ice_clio     : ice   CLIO bulk formulea: compute momentum, heat and freswater fluxes for the sea-ice 
    1718   !!   blk_clio_qsr_oce : shortwave radiation for ocean computed from the cloud cover 
    1819   !!   blk_clio_qsr_ice : shortwave radiation for ice   computed from the cloud cover 
    19    !!   flx_blk_declin : solar declinaison 
     20   !!   flx_blk_declin   : solar declinaison 
    2021   !!---------------------------------------------------------------------- 
    2122   USE oce             ! ocean dynamics and tracers 
     
    4748   INTEGER , PARAMETER ::   jp_vtau = 2           ! index of wind stress (j-component)      (N/m2)    at V-point 
    4849   INTEGER , PARAMETER ::   jp_wndm = 3           ! index of 10m wind module                 (m/s)    at T-point 
    49    INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( - ) 
    50    INTEGER , PARAMETER ::   jp_ccov = 5           ! index of cloud cover                     ( - ) 
     50   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % ) 
     51   INTEGER , PARAMETER ::   jp_ccov = 5           ! index of cloud cover                     ( % ) 
    5152   INTEGER , PARAMETER ::   jp_tair = 6           ! index of 10m air temperature             (Kelvin) 
    5253   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s) 
     
    8182#  include "vectopt_loop_substitute.h90" 
    8283   !!---------------------------------------------------------------------- 
    83    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     84   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    8485   !! $Id$  
    8586   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    9899      !!      the i-component of the stress                (N/m2) 
    99100      !!      the j-component of the stress                (N/m2) 
    100       !!      the 10m wind pseed module                    (m/s) 
     101      !!      the 10m wind speed module                    (m/s) 
    101102      !!      the 10m air temperature                      (Kelvin) 
    102       !!      the 10m specific humidity                    (-) 
    103       !!      the cloud cover                              (-) 
     103      !!      the 10m specific humidity                    (%) 
     104      !!      the cloud cover                              (%) 
    104105      !!      the total precipitation (rain+snow)          (Kg/m2/s) 
    105106      !!              (2) CALL blk_oce_clio 
    106107      !! 
    107108      !!      C A U T I O N : never mask the surface stress fields 
    108       !!                      the stress is assumed to be in the mesh referential 
    109       !!                      i.e. the (i,j) referential 
     109      !!                      the stress is assumed to be in the (i,j) mesh referential 
    110110      !! 
    111111      !! ** Action  :   defined at each time-step at the air-sea interface 
     
    113113      !!              - taum        wind stress module at T-point 
    114114      !!              - wndm        10m wind module at T-point 
    115       !!              - qns, qsr    non-slor and solar heat flux 
    116       !!              - emp, emps   evaporation minus precipitation 
     115      !!              - qns         non-solar heat flux including latent heat of solid  
     116      !!                            precip. melting and emp heat content 
     117      !!              - qsr         solar heat flux 
     118      !!              - emp         upward mass flux (evap. - precip) 
    117119      !!---------------------------------------------------------------------- 
    118120      INTEGER, INTENT( in  ) ::   kt   ! ocean time step 
     
    175177      !                                         ! ====================== ! 
    176178      ! 
    177       CALL fld_read( kt, nn_fsbc, sf )                ! input fields provided at the current time-step 
     179      CALL fld_read( kt, nn_fsbc, sf )                                        ! input fields at the current time-step 
    178180      ! 
    179181#if defined key_lim3       
    180       tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)     !RB ugly patch 
     182      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)                                   ! Tair needed in LIM-3 (!RB ugly patch) 
    181183#endif 
    182       ! 
     184      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_clio( sf, sst_m )      ! surface ocean fluxes using CLIO bulk formulea 
     185      ENDIF                                               !  
     186       
    183187      IF(lwp .AND. nitend-nit000 <= 100 ) THEN 
    184188         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
     
    210214         ENDIF 
    211215      ENDIF 
    212  
    213       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    214           CALL blk_oce_clio( sf, sst_m )                  ! compute the surface ocean fluxes using CLIO bulk formulea 
    215       ENDIF                                               !  
    216216      ! 
    217217   END SUBROUTINE sbc_blk_clio 
     
    239239      !!               - taum        wind stress module at T-point 
    240240      !!               - wndm        10m wind module at T-point 
    241       !!               - qns, qsr    non-slor and solar heat flux 
    242       !!               - emp, emps   evaporation minus precipitation 
     241      !!               - qns         non-solar heat flux including latent heat of solid  
     242      !!                             precip. melting and emp heat content 
     243      !!               - qsr         solar heat flux 
     244      !!               - emp         suface mass flux (evap.-precip.) 
    243245      !!  ** Nota    :   sf has to be a dummy argument for AGRIF on NEC 
    244246      !!---------------------------------------------------------------------- 
     
    257259      REAL(wp) ::   zsst, ztatm, zcco1, zpatm, zcmax, zrmax     !    -         - 
    258260      REAL(wp) ::   zrhoa, zev, zes, zeso, zqatm, zevsqr        !    -         - 
    259       REAL(wp) ::   ztx2, zty2                                  !    -         - 
     261      REAL(wp) ::   ztx2, zty2, zcevap, zcprec                  !    -         - 
    260262      !! 
    261263      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw        ! long-wave heat flux over ocean 
     
    270272      !------------------------------------! 
    271273!CDIR COLLAPSE 
    272       DO jj = 1 , jpj 
    273          DO ji = 1, jpi 
    274             utau(ji,jj) = sf(jp_utau)%fnow(ji,jj) 
    275             vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj) 
    276          END DO 
    277       END DO 
     274      utau(:,:) = sf(jp_utau)%fnow(:,:) 
     275!CDIR COLLAPSE 
     276      vtau(:,:) = sf(jp_vtau)%fnow(:,:) 
     277 
     278      !------------------------------------! 
     279      !   store the wind speed  (wndm )    ! 
     280      !------------------------------------! 
     281!CDIR COLLAPSE 
     282      wndm(:,:) = sf(jp_wndm)%fnow(:,:) 
    278283 
    279284      !------------------------------------! 
     
    291296      CALL lbc_lnk( taum, 'T', 1. ) 
    292297 
    293       !------------------------------------! 
    294       !   store the wind speed  (wndm )    ! 
    295       !------------------------------------! 
    296 !CDIR COLLAPSE 
    297       DO jj = 1 , jpj 
    298          DO ji = 1, jpi 
    299             wndm(ji,jj) = sf(jp_wndm)%fnow(ji,jj) 
    300          END DO 
    301       END DO 
    302  
    303298      !------------------------------------------------! 
    304299      !   Shortwave radiation for ocean and snow/ice   ! 
    305300      !------------------------------------------------! 
    306        
    307301      CALL blk_clio_qsr_oce( qsr ) 
    308302 
     
    401395      !     III    Total FLUXES                                                       ! 
    402396      ! ----------------------------------------------------------------------------- ! 
    403  
    404 !CDIR COLLAPSE 
    405 !CDIR NOVERRCHK 
    406       DO jj = 1, jpj 
    407 !CDIR NOVERRCHK 
    408          DO ji = 1, jpi 
    409             qns (ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj)      ! Downward Non Solar flux 
    410             emp (ji,jj) = zqla(ji,jj) / cevap - sf(jp_prec)%fnow(ji,jj) / rday * tmask(ji,jj,1) 
    411          END DO 
    412       END DO 
    413       emps(:,:) = emp(:,:) 
    414       ! 
     397      zcevap = rcp /  cevap    ! convert zqla ==> evap (Kg/m2/s) ==> m/s ==> W/m2 
     398      zcprec = rcp /  rday     ! convert prec ( mm/day ==> m/s)  ==> W/m2 
     399 
     400!CDIR COLLAPSE 
     401      emp(:,:) = zqla(:,:) / cevap                                        &   ! freshwater flux 
     402         &     - sf(jp_prec)%fnow(:,:) / rday * tmask(:,:,1) 
     403      ! 
     404!CDIR COLLAPSE 
     405      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                        &   ! Downward Non Solar flux 
     406         &     - zqla(:,:)             * pst(:,:)              * zcevap   &   ! remove evap.   heat content at SST in Celcius 
     407         &     + sf(jp_prec)%fnow(:,:) * sf(jp_tair)%fnow(:,:) * zcprec       ! add    precip. heat content at Tair in Celcius 
     408      ! NB: if sea-ice model, the snow precip are computed and the associated heat is added to qns (see blk_ice_clio) 
     409 
    415410      CALL iom_put( "qlw_oce",   zqlw )   ! output downward longwave  heat over the ocean 
    416411      CALL iom_put( "qsb_oce", - zqsb )   ! output downward sensible  heat over the ocean 
     
    425420            &         tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 
    426421      ENDIF 
    427  
     422      ! 
    428423   END SUBROUTINE blk_oce_clio 
    429424 
     
    447442      !! 
    448443      !!  ** Action  :   call albedo_oce/albedo_ice to compute ocean/ice albedo  
    449       !!          computation of snow precipitation 
    450       !!          computation of solar flux at the ocean and ice surfaces 
    451       !!          computation of the long-wave radiation for the ocean and sea/ice 
    452       !!          computation of turbulent heat fluxes over water and ice 
    453       !!          computation of evaporation over water 
    454       !!          computation of total heat fluxes sensitivity over ice (dQ/dT) 
    455       !!          computation of latent heat flux sensitivity over ice (dQla/dT) 
    456       !! 
     444      !!               - snow precipitation 
     445      !!               - solar flux at the ocean and ice surfaces 
     446      !!               - the long-wave radiation for the ocean and sea/ice 
     447      !!               - turbulent heat fluxes over water and ice 
     448      !!               - evaporation over water 
     449      !!               - total heat fluxes sensitivity over ice (dQ/dT) 
     450      !!               - latent heat flux sensitivity over ice (dQla/dT) 
     451      !!               - qns  :  modified the non solar heat flux over the ocean 
     452      !!                         to take into account solid precip latent heat flux 
    457453      !!---------------------------------------------------------------------- 
    458454      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature                   [Kelvin] 
     
    633629      ! 
    634630      ! ----------------------------------------------------------------------------- ! 
    635       !    Total FLUXES                                                       ! 
     631      !    Total FLUXES                                                               ! 
    636632      ! ----------------------------------------------------------------------------- ! 
    637633      ! 
    638634!CDIR COLLAPSE 
    639       p_qns(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - p_qla (:,:,:)      ! Downward Non Solar flux 
     635      p_qns(:,:,:) = z_qlw(:,:,:) - z_qsb(:,:,:) - p_qla(:,:,:)         ! Downward Non Solar flux 
    640636!CDIR COLLAPSE 
    641637      p_tpr(:,:)   = sf(jp_prec)%fnow(:,:) / rday                       ! total precipitation [kg/m2/s] 
     638      ! 
     639      ! ----------------------------------------------------------------------------- ! 
     640      !    Correct the OCEAN non solar flux with the existence of solid precipitation ! 
     641      ! ---------------=====--------------------------------------------------------- ! 
     642!CDIR COLLAPSE 
     643      qns(:,:) = qns(:,:)                                           &   ! update the non-solar heat flux with: 
     644         &     - p_spr(:,:) * lfus                                                  &   ! remove melting solid precip 
     645         &     + p_spr(:,:) * MIN( sf(jp_tair)%fnow(:,:), rt0_snow - rt0 ) * cpic   &   ! add solid P at least below melting 
     646         &     - p_spr(:,:) * sf(jp_tair)%fnow(:,:)                        * rcp        ! remove solid precip. at Tair 
    642647      ! 
    643648!!gm : not necessary as all input data are lbc_lnk... 
     
    667672         CALL prt_ctl(tab2d_1=p_taui , clinfo1=' blk_ice_clio: p_taui : ', tab2d_2=p_tauj , clinfo2=' p_tauj : ') 
    668673      ENDIF 
    669  
    670  
     674      ! 
    671675   END SUBROUTINE blk_ice_clio 
    672676 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r1730 r1859  
    1212   !!            3.0  !  2006-06  (G. Madec) sbc rewritting    
    1313   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put 
     14   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    1415   !!---------------------------------------------------------------------- 
    1516 
     
    4546   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    4647   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point 
    47    INTEGER , PARAMETER ::   jp_humi = 3           ! index of specific humidity               ( - ) 
     48   INTEGER , PARAMETER ::   jp_humi = 3           ! index of specific humidity               ( % ) 
    4849   INTEGER , PARAMETER ::   jp_qsr  = 4           ! index of solar heat                      (W/m2) 
    4950   INTEGER , PARAMETER ::   jp_qlw  = 5           ! index of Long wave                       (W/m2) 
     
    6263   REAL(wp), PARAMETER ::   Cice =    1.63e-3     ! transfer coefficient over ice 
    6364 
    64    !                                !!* Namelist namsbc_core : CORE bulk parameters 
    65    LOGICAL  ::   ln_2m     = .FALSE.     ! logical flag for height of air temp. and hum 
    66    LOGICAL  ::   ln_taudif = .FALSE.     ! logical flag to use the "mean of stress module - module of mean stress" data 
    67    REAL(wp) ::   rn_pfac   = 1.          ! multiplication factor for precipitation 
     65   !                                    !!* Namelist namsbc_core : CORE bulk parameters 
     66   LOGICAL  ::   ln_2m     = .FALSE.     ! air temperature and humidity given at 2m (T) or 10m (F) 
     67   LOGICAL  ::   ln_taudif = .FALSE.     ! (T) use the "mean of stress module - module of mean stress" data or (F) not 
     68   REAL(wp) ::   rn_pfac   = 1.          ! multiplicative factor for precipitation 
    6869 
    6970   !! * Substitutions 
     
    7172#  include "vectopt_loop_substitute.h90" 
    7273   !!---------------------------------------------------------------------- 
    73    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     74   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    7475   !! $Id$ 
    7576   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    8889      !!      the 10m wind velocity (i-component) (m/s)    at T-point 
    8990      !!      the 10m wind velocity (j-component) (m/s)    at T-point 
    90       !!      the specific humidity               ( - ) 
     91      !!      the 10m or 2m specific humidity     ( % ) 
    9192      !!      the solar heat                      (W/m2) 
    9293      !!      the Long wave                       (W/m2) 
    93       !!      the 10m air temperature             (Kelvin) 
     94      !!      the 10m or 2m air temperature       (Kelvin) 
    9495      !!      the total precipitation (rain+snow) (Kg/m2/s) 
    9596      !!      the snow (solid prcipitation)       (kg/m2/s) 
    96       !!   OPTIONAL parameter (see ln_taudif namelist flag): 
    97       !!      the tau diff associated to HF tau   (N/m2)   at T-point  
     97      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T) 
    9898      !!              (2) CALL blk_oce_core 
    9999      !! 
    100100      !!      C A U T I O N : never mask the surface stress fields 
    101       !!                      the stress is assumed to be in the mesh referential 
    102       !!                      i.e. the (i,j) referential 
     101      !!                      the stress is assumed to be in the (i,j) mesh referential 
    103102      !! 
    104103      !! ** Action  :   defined at each time-step at the air-sea interface 
    105104      !!              - utau, vtau  i- and j-component of the wind stress 
    106       !!              - taum        wind stress module at T-point 
    107       !!              - wndm        10m wind module at T-point 
    108       !!              - qns, qsr    non-slor and solar heat flux 
    109       !!              - emp, emps   evaporation minus precipitation 
     105      !!              - taum, wndm  wind stress and 10m wind modules at T-point 
     106      !!              - qns, qsr    non-solar and solar heat flux 
     107      !!              - emp         upward mass flux (evapo. - precip.) 
    110108      !!---------------------------------------------------------------------- 
    111109      INTEGER, INTENT( in  ) ::   kt   ! ocean time step 
    112110      !! 
     111      INTEGER  ::   jf       ! dummy loop indice 
     112      INTEGER  ::   ifld     ! number of files to be read 
    113113      INTEGER  ::   ierror   ! return error code 
    114       INTEGER  ::   ifpr     ! dummy loop indice 
    115       INTEGER  ::   jfld     ! dummy loop arguments 
    116114      !! 
    117115      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files 
    118116      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
    119       TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    120       TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 " 
    121       TYPE(FLD_N) ::   sn_tdif                                 !   "                                 " 
     117      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr             ! informations about the fields to be read 
     118      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow, sn_tdif   !       -                       - 
    122119      NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac,           & 
    123120         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
     
    156153         ! do we use HF tau information? 
    157154         lhftau = ln_taudif 
    158          jfld = jpfld - COUNT( (/.NOT. lhftau/) ) 
     155         ifld = jpfld - COUNT( (/.NOT. lhftau/) ) 
    159156         ! 
    160157         ! set sf structure 
    161          ALLOCATE( sf(jfld), STAT=ierror ) 
     158         ALLOCATE( sf(ifld), STAT=ierror ) 
    162159         IF( ierror > 0 ) THEN 
    163160            CALL ctl_stop( 'sbc_blk_core: unable to allocate sf structure' )   ;   RETURN 
    164161         ENDIF 
    165          DO ifpr= 1, jfld 
    166             ALLOCATE( sf(ifpr)%fnow(jpi,jpj) ) 
    167             ALLOCATE( sf(ifpr)%fdta(jpi,jpj,2) ) 
     162         DO jf = 1, ifld 
     163            ALLOCATE( sf(jf)%fnow(jpi,jpj) ) 
     164            ALLOCATE( sf(jf)%fdta(jpi,jpj,2) ) 
    168165         END DO 
    169166         ! 
     
    173170      ENDIF 
    174171 
     172!!gm    all the below lines should be executed only at nn_fbc frequency, no???   check fldread capability 
     173 
    175174      CALL fld_read( kt, nn_fsbc, sf )                   ! input fields provided at the current time-step 
    176  
     175      ! 
    177176#if defined key_lim3 
    178       tatm_ice(:,:) = sf(jp_tair)%fnow(:,:) 
     177      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)              ! air temperature over ice (LIM3 only) 
    179178#endif 
    180  
    181       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    182           CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m )   ! compute the surface ocean fluxes using CLIO bulk formulea 
    183       ENDIF 
    184       !                                                  ! using CORE bulk formulea 
     179      !                                                  ! surface ocean fluxes using CORE bulk formulea 
     180      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m ) 
     181      ! 
    185182   END SUBROUTINE sbc_blk_core 
    186183    
     
    196193      !!      fields read in sbc_read 
    197194      !!  
    198       !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2) 
    199       !!              - vtau    : j-component of the stress at V-point  (N/m2) 
    200       !!              - taum    : Wind stress module at T-point         (N/m2) 
    201       !!              - wndm    : Wind speed module at T-point          (m/s) 
    202       !!              - qsr     : Solar heat flux over the ocean        (W/m2) 
    203       !!              - qns     : Non Solar heat flux over the ocean    (W/m2) 
    204       !!              - evap    : Evaporation over the ocean            (kg/m2/s) 
    205       !!              - emp(s)  : evaporation minus precipitation       (kg/m2/s) 
     195      !! ** Action  : - utau  : i-component of the stress at U-point  (N/m2) 
     196      !!              - vtau  : j-component of the stress at V-point  (N/m2) 
     197      !!              - taum  : Wind stress module at T-point         (N/m2) 
     198      !!              - wndm  : 10m Wind speed module at T-point      (m/s) 
     199      !!              - qsr   : Solar heat flux over the ocean        (W/m2) 
     200      !!              - qns   : Non Solar heat flux over the ocean    (W/m2) 
     201      !!                        including the latent heat of solid  
     202      !!                        precip. melting and emp heat content 
     203      !!              - emp   : upward mass flux (evap. - precip.)    (kg/m2/s) 
    206204      !! 
    207205      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC 
    208206      !!--------------------------------------------------------------------- 
    209       TYPE(fld), INTENT(in), DIMENSION(:)       ::   sf    ! input data 
    210       REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pst   ! surface temperature                      [Celcius] 
    211       REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pu    ! surface current at U-point (i-component) [m/s] 
    212       REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pv    ! surface current at V-point (j-component) [m/s] 
    213  
     207      TYPE(fld), INTENT(in), DIMENSION(:)       ::   sf    ! input data (forcing field structure) 
     208      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pst   ! surface temperature                      [Celcius] 
     209      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pu    ! surface current at U-point (i-component) [m/s] 
     210      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pv    ! surface current at V-point (j-component) [m/s] 
     211      !! 
    214212      INTEGER  ::   ji, jj     ! dummy loop indices 
    215       REAL(wp) ::   zcoef_qsatw 
    216       REAL(wp) ::   zztmp                                 ! temporary variable 
     213      REAL(wp) ::   zcoef_qsatw, zztmp                    ! temporary scalar 
    217214      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    218215      REAL(wp), DIMENSION(jpi,jpj) ::   zqsatw            ! specific humidity at pst 
     
    230227      zcoef_qsatw = 0.98 * 640380. / rhoa 
    231228       
    232       zst(:,:) = pst(:,:) + rt0      ! converte Celcius to Kelvin (and set minimum value far above 0 K) 
     229      zst(:,:) = pst(:,:) + rt0      ! converte SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    233230 
    234231      ! ----------------------------------------------------------------------------- ! 
     
    262259      ! ocean albedo assumed to be 0.066 
    263260!CDIR COLLAPSE 
    264       qsr (:,:) = ( 1. - 0.066 ) * sf(jp_qsr)%fnow(:,:) * tmask(:,:,1)                                 ! Short Wave 
     261      qsr (:,:) = ( 1. - 0.066 ) * sf(jp_qsr)%fnow(:,:) * tmask(:,:,1)                                     ! Short Wave 
    265262!CDIR COLLAPSE 
    266263      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
     
    353350      
    354351!CDIR COLLAPSE 
    355       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)      ! Downward Non Solar flux 
    356 !CDIR COLLAPSE 
    357       emp (:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:) * rn_pfac * tmask(:,:,1) 
    358 !CDIR COLLAPSE 
    359       emps(:,:) = emp(:,:) 
     352      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
     353         &         - sf(jp_prec)%fnow(:,:) * rn_pfac  ) * tmask(:,:,1) 
     354!CDIR COLLAPSE 
     355      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                       &   ! Downward Non Solar flux 
     356         &     - sf(jp_snow)%fnow(:,:) * lfus                            &   ! remove latent melting heat for solid precip 
     357         &     - zevap(:,:) * pst(ji,jj) * rcp                           &   ! remove evap heat content at SST 
     358         &     + ( sf(jp_prec)%fnow(:,:) - sf(jp_snow)%fnow(:,:) )       &   ! add liquid precip heat content at Tair 
     359         &     * ( sf(jp_tair)%fnow(:,:) - rt0 ) * rcp                   &    
     360         &     + sf(jp_snow)%fnow(:,:)                                   &   ! add solid  precip heat content at min(Tair,Tsnow) 
     361         &     * ( MIN( sf(jp_tair)%fnow(:,:), rt0_snow ) - rt0 ) * cpic  
    360362      ! 
    361363      CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean 
     
    392394      !! caution : the net upward water flux has with mm/day unit 
    393395      !!--------------------------------------------------------------------- 
    394       REAL(wp), INTENT(in   ), DIMENSION(:,:,:)      ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin] 
    395       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)    ::   pui      ! ice surface velocity (i- and i- components      [m/s] 
    396       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)    ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid) 
    397       REAL(wp), INTENT(in   ), DIMENSION(:,:,:)      ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%] 
    398       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_taui   ! i- & j-components of surface ice stress        [N/m2] 
    399       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid) 
    400       REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2] 
    401       REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2] 
    402       REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2] 
    403       REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2] 
    404       REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2] 
    405       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s] 
    406       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s] 
    407       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%] 
    408       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%] 
    409       CHARACTER(len=1), INTENT(in   )                ::   cd_grid  ! ice grid ( C or B-grid) 
    410       INTEGER, INTENT(in   )                         ::   pdim     ! number of ice categories 
     396      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin] 
     397      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pui      ! ice surface velocity (i- and i- components      [m/s] 
     398      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid) 
     399      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%] 
     400      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2] 
     401      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid) 
     402      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2] 
     403      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2] 
     404      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2] 
     405      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2] 
     406      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2] 
     407      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s] 
     408      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s] 
     409      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%] 
     410      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%] 
     411      CHARACTER(len=1), INTENT(in   )             ::   cd_grid  ! ice grid ( C or B-grid) 
     412      INTEGER, INTENT(in   )                      ::   pdim     ! number of ice categories 
    411413      !! 
    412414      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    413415      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
    414       REAL(wp) ::   zst2, zst3 
    415       REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb 
    416       REAL(wp) ::   zcoef_frca                                   ! fractional cloud amount 
    417       REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                  ! relative wind module and components at F-point 
    418       REAL(wp) ::             zwndi_t , zwndj_t                  ! relative wind components at T-point 
    419       REAL(wp), DIMENSION(jpi,jpj) ::   z_wnds_t                 ! wind speed ( = | U10m - U_ice | ) at T-point 
    420       REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_qlw               ! long wave heat flux over ice 
    421       REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_qsb               ! sensible  heat flux over ice 
    422       REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_dqlw              ! long wave heat sensitivity over ice 
    423       REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_dqsb              ! sensible  heat sensitivity over ice 
     416      REAL(wp) ::   zst2, zcoef_wnorm , zcoef_dqlw              ! 
     417      REAL(wp) ::   zst3, zcoef_wnorm2, zcoef_dqla, zcoef_dqsb  ! 
     418      REAL(wp) ::   zcoef_frca                                  ! fractional cloud amount 
     419      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                 ! relative wind module and components at F-point 
     420      REAL(wp) ::             zwndi_t , zwndj_t                 ! relative wind components at T-point 
     421      REAL(wp), DIMENSION(jpi,jpj)      ::   z_wnds_t           ! wind speed ( = | U10m - U_ice | ) at T-point 
     422      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_qlw              ! long wave heat flux over ice 
     423      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_qsb              ! sensible  heat flux over ice 
     424      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_dqlw             ! long wave heat sensitivity over ice 
     425      REAL(wp), DIMENSION(jpi,jpj,pdim) ::   z_dqsb             ! sensible  heat sensitivity over ice 
    424426      !!--------------------------------------------------------------------- 
    425427 
     
    576578         CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ') 
    577579      ENDIF 
    578  
     580      ! 
    579581   END SUBROUTINE blk_ice_core 
    580582   
    581583 
    582584   SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a,   & 
    583       &                        dU, Cd, Ch, Ce   ) 
     585      &                        dU , Cd , Ch   , Ce   ) 
    584586      !!---------------------------------------------------------------------- 
    585587      !!                      ***  ROUTINE  turb_core  *** 
     
    704706 
    705707 
    706     SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu) 
     708    SUBROUTINE TURB_CORE_2Z( zt  , zu, sst, T_zt, q_sat,   & 
     709      &                      q_zt, dU, Cd , Ch  , Ce   , T_zu, q_zu) 
    707710      !!---------------------------------------------------------------------- 
    708711      !!                      ***  ROUTINE  turb_core  *** 
     
    838841         Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct 
    839842         !! 
    840          !! 
    841843      END DO 
    842       !! 
     844      ! 
    843845    END SUBROUTINE TURB_CORE_2Z 
    844846 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbccpl.F90

    r1833 r1859  
    44   !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode 
    55   !!====================================================================== 
    6    !! History :  2.0  !  06-2007  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod 
    7    !!            3.0  !  02-2008  (G. Madec, C Talandier)  surface module 
    8    !!            3.1  !  02-2009  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface 
     6   !! History :  2.0  !  2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod 
     7   !!            3.0  !  2008-02  (G. Madec, C Talandier)  surface module 
     8   !!            3.1  !  2009-02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface 
     9   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_oasis3 || defined key_oasis4 
     
    156157#  include "vectopt_loop_substitute.h90" 
    157158   !!---------------------------------------------------------------------- 
    158    !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     159   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    159160   !! $Id$ 
    160161   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    562563      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid  
    563564      !!                        taum, wndm   wind stres and wind speed module at T-point 
    564       !!                        qns , qsr    non solar and solar ocean heat fluxes   ('ocean only case) 
    565       !!                        emp = emps   evap. - precip. (- runoffs) (- calving) ('ocean only case) 
     565      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case) 
     566      !!                                     and the latent heat flux of solid precip. melting 
     567      !!                        qsr          solar ocean heat fluxes   (ocean only case) 
     568      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 
    566569      !!---------------------------------------------------------------------- 
    567570      INTEGER, INTENT(in) ::   kt       ! ocean model time step index 
     
    697700      ENDIF 
    698701 
    699       ! u(v)tau and taum will be modified by ice model (wndm will be changed by PISCES) 
     702      ! u(v)tau and taum will be modified by ice model (and wndm will be changed by PISCES) 
    700703      ! -> need to be reset before each call of the ice/fsbc       
    701704      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN 
     
    712715         !                                                   ! ========================= ! 
    713716         ! 
    714          !                                                       ! non solar heat flux over the ocean (qns) 
    715          IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(:,:,jpr_qnsoce) 
    716          IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(:,:,jpr_qnsmix)         
    717          !   energy for melting solid precipitation over free ocean 
    718          zcoef = xlsn / rhosn 
    719          qns(:,:) = qns(:,:) - frcv(:,:,jpr_snow) * zcoef 
    720          !                                                       ! solar flux over the ocean          (qsr) 
    721          IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(:,:,jpr_qsroce)  
    722          IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(:,:,jpr_qsrmix) 
    723          ! 
    724          !                                                       ! total freshwater fluxes over the ocean (emp, emps) 
     717         !                                                       ! total freshwater fluxes over the ocean (emp) 
    725718         SELECT CASE( TRIM( cn_rcv_emp ) )                                    ! evaporation - precipitation 
    726719         CASE( 'conservative' ) 
     
    752745!!         ENDIF 
    753746!!gm  end of internal cooking 
    754          ! 
    755          emps(:,:) = emp(:,:)                                        ! concentration/dilution = emp 
    756    
    757          !                                                           ! 10 m wind speed 
    758          IF( srcv(jpr_w10m)%laction )   wndm(:,:) = frcv(:,:,jpr_w10m) 
     747         !   
     748         !                                                       ! non solar heat flux over the ocean (qns) 
     749         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(:,:,jpr_qnsoce) 
     750         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(:,:,jpr_qnsmix) 
     751         ! 
     752         zcoef = xlsn / rhosn                                    ! qns update over free ocean with: 
     753         qns(:,:) = qns(:,:) - frcv(:,:,jpr_snow) * zcoef            ! energy for melting solid precipitation over free ocean 
     754            &                - emp(:,:) * sst_m(:,:) * rcp           ! remove heat content due to mass flux (assumed to be at SST) 
     755         ! 
     756         !                                                       ! solar flux over the ocean          (qsr) 
     757         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(:,:,jpr_qsroce)  
     758         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(:,:,jpr_qsrmix) 
    759759         ! 
    760760#if defined  key_cpl_carbon_cycle 
    761          !                                                              ! atmosph. CO2 (ppm) 
     761         !                                                       ! atmosph. CO2 (ppm) 
    762762         IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(:,:,jpr_co2) 
    763763#endif 
    764  
     764         ! 
    765765      ENDIF 
    766766      ! 
     
    10461046      !!---------------------------------------------------------------------- 
    10471047      zicefr(:,:,1) = 1.- p_frld(:,:,1) 
    1048       IF( lk_diaar5 )   zcptn(:,:) = rcp * tn(:,:,1) 
     1048      zcptn(:,:) = rcp * sst_m(:,:) 
    10491049      ! 
    10501050      !                                                      ! ========================= ! 
     
    11181118            &                                                   +          pist(:,:,1)   * zicefr(:,:,1) ) ) 
    11191119      END SELECT 
    1120       !                                                           ! snow melting heat flux .... 
    1121       !   energy for melting solid precipitation over ice-free ocean 
    1122       zcoef = xlsn / rhosn 
     1120      ! 
     1121      zcoef = xlsn / rhosn                                        ! qns_tot update over free ocean with: 
    11231122      ztmp(:,:) = p_frld(:,:,1) * zsnow(:,:) * zcoef 
    1124       pqns_tot(:,:) = pqns_tot(:,:) - ztmp(:,:) 
    1125       IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + zsnow(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average) 
    1126 !!gm 
    1127 !!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in  
    1128 !!    the flux that enter the ocean.... 
    1129 !!    moreover 1 - it is not diagnose anywhere....  
    1130 !!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not... 
    1131 !! 
    1132 !! similar job should be done for snow and precipitation temperature 
    1133       !                                                           ! Iceberg melting heat flux .... 
    1134       !   energy for iceberg melting  
    1135       IF( srcv(jpr_cal)%laction ) THEN  
     1123      pqns_tot(:,:) = pqns_tot(:,:)                       & 
     1124         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting 
     1125         &          + (  pemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST) 
     1126         &             - pemp_ice(:,:) * p_frld(:,:,1)  ) * zcptn(:,:)  
     1127         ! 
     1128      IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp(:,:) + zsnow(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average) 
     1129 
     1130!!gm BUG ???   just above the right value should be : ztmp + zsnow*p_frld*zcptn 
     1131 
     1132      !                                
     1133      IF( srcv(jpr_cal)%laction ) THEN                                 ! remove the latent heat flux of iceberg melting 
    11361134         zcoef = xlic / rhoic 
    11371135         ztmp(:,:) = frcv(:,:,jpr_cal) * zcoef 
    11381136         pqns_tot(:,:) = pqns_tot(:,:) - ztmp(:,:) 
     1137         ! 
    11391138         IF( lk_diaar5 )   CALL iom_put( 'hflx_cal_cea', ztmp + frcv(:,:,jpr_cal) * zcptn(:,:) )   ! heat flux from calving 
    11401139      ENDIF 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcflx.F90

    r1730 r1859  
    44   !! Ocean forcing:  momentum, heat and freshwater flux formulation 
    55   !!===================================================================== 
    6    !! History :  9.0   !  06-06  (G. Madec)  Original code 
     6   !! History :  1.0  !  2006-06  (G. Madec)  Original code 
     7   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    78   !!---------------------------------------------------------------------- 
    89 
     
    5253#  include "vectopt_loop_substitute.h90" 
    5354   !!---------------------------------------------------------------------- 
    54    !!   OPA 9.0 , LOCEAN-IPSL (2006)  
     55   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    5556   !! $Id$ 
    5657   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    7475      !! 
    7576      !!      CAUTION :  - never mask the surface stress fields 
    76       !!                 - the stress is assumed to be in the mesh referential 
    77       !!                   i.e. the (i,j) referential 
     77      !!                 - the stress is assumed to be in the (i,j) mesh referential 
    7878      !! 
    7979      !! ** Action  :   update at each time-step 
     
    8181      !!              - taum        wind stress module at T-point 
    8282      !!              - wndm        10m wind module at T-point 
    83       !!              - qns, qsr    non-slor and solar heat flux 
    84       !!              - emp, emps   evaporation minus precipitation 
     83      !!              - qns         non solar heat flux including heat flux due to emp 
     84      !!              - qsr         solar heat flux 
     85      !!              - emp         upward mass flux (evap. - precip.) 
    8586      !!---------------------------------------------------------------------- 
    8687      INTEGER, INTENT(in) ::   kt   ! ocean time step 
     
    136137      ENDIF 
    137138 
    138       CALL fld_read( kt, nn_fsbc, sf )           ! Read input fields and provides the 
    139       !                                          ! input fields at the current time-step 
     139      CALL fld_read( kt, nn_fsbc, sf )           ! input fields at the current time-step 
    140140 
    141       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    142          ! 
    143          ! set the ocean fluxes from read fields 
     141      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN       ! set the ocean fluxes from read fields 
    144142!CDIR COLLAPSE 
    145143         DO jj = 1, jpj 
     
    151149               emp (ji,jj) = sf(jp_emp )%fnow(ji,jj) 
    152150            END DO 
    153          END DO 
    154           
    155          ! module of wind stress and wind speed at T-point 
    156          zcoef = 1. / ( zrhoa * zcdrag )  
     151         END DO          
     152         !                                       ! add to qns the heat due to e-p 
     153         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp      ! mass flux are at SST 
     154 
     155         zcoef = 1. / ( zrhoa * zcdrag )         ! module of wind stress and wind speed at T-point 
    157156!CDIR NOVERRCHK 
    158157         DO jj = 2, jpjm1 
     
    167166         END DO 
    168167         CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
    169  
    170          ! Initialization of emps (when no ice model) 
    171          emps(:,:) = emp (:,:)  
    172168                   
    173          ! control print (if less than 100 time-step asked) 
    174          IF( nitend-nit000 <= 100 .AND. lwp ) THEN 
     169         IF( nitend-nit000 <= 100 .AND. lwp ) THEN      ! control print (if less than 100 time-step asked) 
    175170            WRITE(numout,*)  
    176171            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK' 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r1822 r1859  
    44   !! Ocean fluxes   : domain averaged freshwater budget 
    55   !!====================================================================== 
    6    !! History :  8.2  !  01-02  (E. Durand)  Original code 
    7    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
    8    !!            9.0  !  06-08  (G. Madec)  Surface module 
    9    !!            9.2  !  09-07  (C. Talandier) emp mean s spread over erp area  
     6   !! History :  OPA  !  2001-02  (E. Durand)  Original code 
     7   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
     8   !!            3.0  !  2006-08  (G. Madec)  Surface module 
     9   !!            3.2  !  2009-07  (C. Talandier) emp mean s spread over erp area  
     10   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) embedded sea-ice case 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    2728   PRIVATE 
    2829 
    29    PUBLIC   sbc_fwb      ! routine called by step 
    30  
    31    REAL(wp) ::   a_fwb_b            ! annual domain averaged freshwater budget 
    32    REAL(wp) ::   a_fwb              ! for 2 year before (_b) and before year. 
    33    REAL(wp) ::   empold             ! empold to be suppressed 
    34    REAL(wp) ::   area               ! global mean ocean surface (interior domain) 
     30   PUBLIC   sbc_fwb    ! routine called by step 
     31 
     32   REAL(wp) ::   a_fwb_b   ! annual domain averaged freshwater budget 
     33   REAL(wp) ::   a_fwb     ! for 2 year before (_b) and before year. 
     34   REAL(wp) ::   empold    ! empold to be suppressed 
     35   REAL(wp) ::   area      ! global mean ocean surface (interior domain) 
    3536 
    3637   REAL(wp), DIMENSION(jpi,jpj) ::   e1e2_i    ! area of the interior domain (e1t*e2t*tmask_i) 
     
    4041#  include "vectopt_loop_substitute.h90" 
    4142   !!---------------------------------------------------------------------- 
    42    !!  OPA 9.0 , LOCEAN-IPSL (2006)  
     43   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    4344   !! $Id$ 
    4445   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    5758      !!                =2 annual global mean corrected from previous year 
    5859      !!                =3 global mean of emp set to zero at each nn_fsbc time step 
    59       !!                   & spread out over erp area depending its sign 
     60      !!                   and spread out over erp area depending its sign 
    6061      !!---------------------------------------------------------------------- 
    6162      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index 
     
    6364      INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index 
    6465      !! 
    65       INTEGER  ::   inum                  ! temporary logical unit 
    66       INTEGER  ::   ikty, iyear           !  
    67       REAL(wp) ::   z_emp, z_emp_nsrf, zsum_emp, zsum_erp       ! temporary scalars 
    68       REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread 
    69       REAL(wp), DIMENSION(jpi,jpj) ::   ztmsk_neg, ztmsk_pos, ztmsk_tospread 
    70       REAL(wp), DIMENSION(jpi,jpj) ::   z_wgt, zerp_cor 
     66      INTEGER  ::   inum          ! temporary logical unit 
     67      INTEGER  ::   ikty, iyear   !  
     68      REAL(wp) ::   z_emp, z_emp_nsrf, zsum_emp, zsum_erp, zcoef     ! temporary scalars 
     69      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread             !    -         - 
     70      REAL(wp), DIMENSION(jpi,jpj) ::   ztmsk_neg, ztmsk_tospread    ! 2D workspace 
     71      REAL(wp), DIMENSION(jpi,jpj) ::   ztmsk_pos, z_wgt, zerp_cor   !  -      - 
    7172      !!---------------------------------------------------------------------- 
    7273      ! 
     
    8384            IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   & 
    8485               &   CALL ctl_stop( 'The option nn_fwb = 3 must be associated to nn_sssr = 2 ' ) 
    85              
    8686         ENDIF 
    8787         ! 
     
    9999         CALL ctl_stop( ctmp1 ) 
    100100         ! 
    101           
    102       ! 
    103101      CASE ( 1 )                               ! global mean emp set to zero 
    104102         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 
    105103            z_emp = SUM( e1e2_i(:,:) * emp(:,:) ) / area 
    106104            IF( lk_mpp )   CALL  mpp_sum( z_emp    )   ! sum over the global domain 
    107             emp (:,:) = emp (:,:) - z_emp 
    108             emps(:,:) = emps(:,:) - z_emp 
     105            zcoef = z_emp * rcp 
     106            emp(:,:) = emp(:,:) - z_emp 
     107            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) 
    109108         ENDIF 
    110109         ! 
     
    138137         ! correct the freshwater fluxes 
    139138         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 
    140             emp (:,:) = emp (:,:) + empold 
    141             emps(:,:) = emps(:,:) + empold 
     139            zcoef = z_emp * rcp 
     140            emp(:,:) = emp (:,:) + empold 
     141            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) 
    142142         ENDIF 
    143143         ! 
     
    152152         ! 
    153153         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 
    154             ! Select <0 and >0 area of erp 
     154            !                                       ! Select <0 and >0 area of erp 
    155155            ztmsk_pos(:,:) = tmask_i(:,:) 
    156             WHERE( erp < 0.e0 ) ztmsk_pos = 0.e0 
     156            WHERE( erp(:,:) < 0.e0 )   ztmsk_pos(:,:) = 0.e0 
    157157            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:) 
    158158 
    159             ! Area filled by <0 and >0 erp  
    160             zsurf_neg = SUM( e1e2_i(:,:)*ztmsk_neg(:,:) ) 
    161             zsurf_pos = SUM( e1e2_i(:,:)*ztmsk_pos(:,:) ) 
     159            !                                       ! Area filled by <0 and >0 erp  
     160            zsurf_neg = SUM(  e1e2_i(:,:) * ztmsk_neg(:,:) ) 
     161            zsurf_pos = SUM(  e1e2_i(:,:) * ztmsk_pos(:,:) ) 
    162162         
    163             ! emp global mean  
     163            !                                       ! emp global mean  
    164164            z_emp = SUM( e1e2_i(:,:) * emp(:,:) ) / area 
    165165            ! 
    166             IF( lk_mpp )   CALL  mpp_sum( z_emp ) 
     166            IF( lk_mpp )   CALL  mpp_sum( z_emp     ) 
    167167            IF( lk_mpp )   CALL  mpp_sum( zsurf_neg ) 
    168168            IF( lk_mpp )   CALL  mpp_sum( zsurf_pos ) 
    169169             
    170             IF( z_emp < 0.e0 ) THEN 
    171                 ! to spread out over >0 erp area to increase evaporation damping process 
     170            IF( z_emp < 0.e0 ) THEN                 ! spread out over >0 erp area to increase evaporation damping process 
    172171                zsurf_tospread = zsurf_pos 
    173172                ztmsk_tospread(:,:) = ztmsk_pos(:,:) 
    174             ELSE 
    175                 ! to spread out over <0 erp area to increase precipitation damping process 
     173            ELSE                                    ! spread out over <0 erp area to increase precipitation damping process 
    176174                zsurf_tospread = zsurf_neg 
    177175                ztmsk_tospread(:,:) = ztmsk_neg(:,:) 
     
    192190            CALL lbc_lnk( zerp_cor, 'T', 1. ) 
    193191 
    194             emp (:,:) = emp (:,:) + zerp_cor(:,:) 
    195             emps(:,:) = emps(:,:) + zerp_cor(:,:) 
    196             erp (:,:) = erp (:,:) + zerp_cor(:,:) 
     192            emp(:,:) = emp(:,:) + zerp_cor(:,:) 
     193            qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:) 
     194            erp(:,:) = erp(:,:) + zerp_cor(:,:) 
    197195             
    198196            IF( nprint == 1 .AND. lwp ) THEN 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcice_if.F90

    r1730 r1859  
    3030#  include "domzgr_substitute.h90" 
    3131   !!---------------------------------------------------------------------- 
    32    !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     32   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    3333   !! $Id$ 
    3434   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    7272         !             !   file    ! frequency !  variable  ! time intep !  clim  ! 'yearly' or ! weights  ! rotation   ! 
    7373         !             !   name    !  (hours)  !   name     !   (T/F)    !  (T/F) !  'monthly'  ! filename ! pairs      !  
    74          sn_ice = FLD_N('ice_cover',    -1    ,  'ice_cov' ,  .true.    , .true. ,   'yearly'  , ''       , ''         ) 
    75  
    76          REWIND ( numnam )               ! ... read in namlist namiif 
     74         sn_ice = FLD_N('ice_cover',    -1     ,  'ice_cov' ,  .true.    , .true. ,   'yearly'  , ''       , ''         ) 
     75         ! 
     76         REWIND ( numnam )               ! read in namlist namiif 
    7777         READ   ( numnam, namsbc_iif ) 
    78  
     78         ! 
    7979         ALLOCATE( sf_ice(1), STAT=ierror ) 
    8080         IF( ierror > 0 ) THEN 
     
    8383         ALLOCATE( sf_ice(1)%fnow(jpi,jpj) ) 
    8484         ALLOCATE( sf_ice(1)%fdta(jpi,jpj,2) ) 
    85  
    86  
    87          ! fill sf_ice with sn_ice and control print 
     85         ! 
     86         !                               ! fill sf_ice with sn_ice and control print 
    8887         CALL fld_fill( sf_ice, (/ sn_ice /), cn_dir, 'sbc_ice_if', 'ice-if sea-ice model', 'namsbc_iif' ) 
    8988         ! 
    9089      ENDIF 
    9190 
    92       CALL fld_read( kt, nn_fsbc, sf_ice )           ! Read input fields and provides the 
    93       !                                              ! input fields at the current time-step 
     91      CALL fld_read( kt, nn_fsbc, sf_ice )           ! Read input fields at the current time-step 
    9492       
    9593      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 
     
    10199         fr_i(:,:) = tfreez( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius] 
    102100 
    103          ! Flux and ice fraction computation 
    104101!CDIR COLLAPSE 
    105          DO jj = 1, jpj 
     102         DO jj = 1, jpj          ! Flux and ice fraction computation 
    106103            DO ji = 1, jpi 
    107104               ! 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcmod.F90

    r1792 r1859  
    44   !! Surface module :  provide to the ocean its surface boundary condition 
    55   !!====================================================================== 
    6    !! History :  3.0   !  07-2006  (G. Madec)  Original code 
    7    !!             -    !  08-2008  (S. Masson, E. .... ) coupled interface 
     6   !! History :  3.0  !  2007-06  (G. Madec)  Original code 
     7   !!            3.1  !  2008-08  (S. Masson, E. .... ) coupled interface 
     8   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    89   !!---------------------------------------------------------------------- 
    910 
     
    4950#  include "domzgr_substitute.h90" 
    5051   !!---------------------------------------------------------------------- 
    51    !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     52   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    5253   !! $Id$ 
    5354   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    8687!!gm  IF( lk_sbc_cpl       ) THEN   ;   ln_cpl      = .TRUE.   ;   ELSE   ;   ln_cpl      = .FALSE.   ;   ENDIF 
    8788 
    88       IF ( Agrif_Root() ) THEN 
     89      IF( Agrif_Root() ) THEN 
    8990        IF( lk_lim2 )            nn_ice      = 2 
    9091        IF( lk_lim3 )            nn_ice      = 3 
     
    123124      ENDIF 
    124125      IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero 
     126 
     127      emps(:,:) = 0.e0                             ! the salt flux will be computed (i.e. will be non-zero) only if  
     128      !                                            ! sea-ice is present, or lk_vvl=F, or surface salt restoring is used. 
    125129 
    126130      !                                            ! restartability    
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r1730 r1859  
    88   !!            3.0  !  2006-07  (G. Madec)  Surface module  
    99   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put 
     10   !!             -   !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    4142   REAL(wp), PUBLIC, DIMENSION(jpk)     ::   rnfmsk_z    !: river mouth mask (vert.) 
    4243 
     44   REAL(wp) ::   rfact_rcp   ! = rn_rfact * rcp 
    4345   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnf   ! structure of input SST (file information, fields read) 
    4446 
    4547   !!---------------------------------------------------------------------- 
    46    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     48   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    4749   !! $Id$ 
    4850   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    6466      !!---------------------------------------------------------------------- 
    6567      INTEGER, INTENT(in) ::   kt          ! ocean time step 
    66       !! 
     68      ! 
    6769      INTEGER  ::   ji, jj   ! dummy loop indices 
    6870      INTEGER  ::   ierror   ! temporary integer 
     
    7880            ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,2) ) 
    7981         ENDIF 
    80          CALL sbc_rnf_init(sf_rnf) 
     82         CALL sbc_rnf_init( sf_rnf ) 
     83         ! 
     84         rfact_rcp = rn_rfact * rcp 
    8185      ENDIF 
    8286 
     
    8589         !                                                !-------------------! 
    8690         ! 
    87          CALL fld_read( kt, nn_fsbc, sf_rnf )   ! Read Runoffs data and provides it 
    88          !                                      ! at the current time-step 
    89  
    90          ! Runoff reduction only associated to the ORCA2_LIM configuration 
    91          ! when reading the NetCDF file runoff_1m_nomask.nc 
    92          IF( cp_cfg == 'orca' .AND. jp_cfg == 2 )   THEN 
    93             DO jj = 1, jpj 
     91         CALL fld_read( kt, nn_fsbc, sf_rnf )                    ! Read Runoffs data at the current time-step 
     92         ! 
     93!!gm CAUTION this is ugly  ===>>> to be removed! 
     94         IF( cp_cfg == 'orca' .AND. jp_cfg == 2 )   THEN         ! Runoff reduction only associated to the ORCA2_LIM configuration 
     95            DO jj = 1, jpj                                       ! when reading the NetCDF file runoff_1m_nomask.nc 
    9496               DO ji = 1, jpi 
    9597                  IF( gphit(ji,jj) > 40 .AND. gphit(ji,jj) < 65 )   sf_rnf(1)%fnow(ji,jj) = 0.85 * sf_rnf(1)%fnow(ji,jj) 
     
    9799            END DO 
    98100         ENDIF 
    99  
    100          ! C a u t i o n : runoff is negative and in kg/m2/s  
    101  
    102          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    103             emp (:,:) = emp (:,:) - rn_rfact * ABS( sf_rnf(1)%fnow(:,:) ) 
    104             emps(:,:) = emps(:,:) - rn_rfact * ABS( sf_rnf(1)%fnow(:,:) ) 
     101         ! 
     102         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                    ! C a u t i o n : runoff is negative and in kg/m2/s  
     103            emp(:,:) = emp(:,:) - rn_rfact  * ABS( sf_rnf(1)%fnow(:,:) )                ! mass flux 
     104            qns(:,:) = qns(:,:) + rfact_rcp * ABS( sf_rnf(1)%fnow(:,:) ) * sst_m(:,:)   ! its associated heat content (at SST) 
     105            ! 
    105106            CALL iom_put( "runoffs", sf_rnf(1)%fnow )         ! runoffs 
    106107         ENDIF 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbcssr.F90

    r1730 r1859  
    66   !! History :  3.0  !  2006-06  (G. Madec)  Original code 
    77   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put 
     8   !!             -   !  2009-07  (C. Talandier, G. Madec)  Add a bound to the Erp  
     9   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    810   !!---------------------------------------------------------------------- 
    911 
     
    3739   REAL(wp)        ::   rn_deds     = -27.70    ! restoring factor on SST and SSS 
    3840   LOGICAL         ::   ln_sssr_bnd = .false.   ! flag to bound erp term  
    39    REAL(wp)        ::   rn_sssr_bnd =   0.e0    ! ABS(Max./Min.) value of erp term [mm/day] 
     41   REAL(wp)        ::   rn_sssr_bnd =   4.e0    ! ABS(Max./Min.) value of erp term [mm/day] 
    4042 
    4143   REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange 
     
    4648#  include "domzgr_substitute.h90" 
    4749   !!---------------------------------------------------------------------- 
    48    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     50   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    4951   !! $Id$ 
    5052   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    7072      !! 
    7173      INTEGER  ::   ji, jj   ! dummy loop indices 
    72       REAL(wp) ::   zerp     ! local scalar for evaporation damping 
    73       REAL(wp) ::   zqrp     ! local scalar for heat flux damping 
    74       REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor 
    75       REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor 
    7674      INTEGER  ::   ierror   ! return error code 
    77       !! 
    78       CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files 
     75      REAL(wp) ::   zerp, zqrp, zsrp, zerp_bnd    ! local scalar 
     76      !! 
     77      CHARACTER(len=100) ::  cn_dir = './'   ! Root directory for location of ssr files 
    7978      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read 
    8079      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd 
     
    8584         !                                            ! -------------------- ! 
    8685         !                            !* set file information 
    87          cn_dir  = './'            ! directory in which the model is executed 
    8886         ! ... default values (NB: frequency positive => hours, negative => months) 
    8987         !            !   file    ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! 
     
    158156                  END DO 
    159157               END DO 
    160                CALL iom_put( "qrp", qrp )                             ! heat flux damping 
    161             ENDIF 
    162             ! 
    163             IF( nn_sssr == 1 ) THEN                   !* Salinity damping term (salt flux, emps only) 
     158               CALL iom_put( "qrp", qrp )                             ! heat flux damping  
     159            ENDIF 
     160            ! 
     161            IF( nn_sssr == 1 ) THEN                   !* Salinity damping term (salt flux only (emps)) 
    164162               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s] 
    165163!CDIR COLLAPSE 
     
    167165                  DO ji = 1, jpi 
    168166                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    169                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj) )   & 
    170                         &        / ( sss_m(ji,jj) + 1.e-20   ) 
     167                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj) ) 
    171168                     emps(ji,jj) = emps(ji,jj) + zerp 
    172                      erp( ji,jj) = zerp 
     169                     erp( ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20  )  ! converted into an equivalent emp (diag. only) 
    173170                  END DO 
    174171               END DO 
    175172               CALL iom_put( "erp", erp )                             ! freshwater flux damping 
    176173               ! 
    177             ELSEIF( nn_sssr == 2 ) THEN               !* Salinity damping term (volume flux, emp and emps) 
     174            ELSEIF( nn_sssr == 2 ) THEN               !* Salinity damping term (volume flux (emp) and qns) 
    178175               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s] 
    179176               zerp_bnd = rn_sssr_bnd / rday                          !       -              -     
     
    183180                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    184181                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj) )   & 
    185                         &        / ( sss_m(ji,jj) + 1.e-20   ) 
     182                        &        / MAX(  sss_m(ji,jj), 1.e-20  ) 
    186183                     IF( ln_sssr_bnd )   zerp = SIGN( 1., zerp ) * MIN( zerp_bnd, ABS(zerp) ) 
    187                      emp (ji,jj) = emp (ji,jj) + zerp 
    188                      emps(ji,jj) = emps(ji,jj) + zerp 
    189                      erp (ji,jj) = zerp 
     184!!gm better coding   IF( ln_sssr_bnd )   zerp = MAX( -zerp_bnd, MIN( zerp, zerp_bnd )  ) 
     185                     emp(ji,jj) = emp(ji,jj) + zerp 
     186                     qns(ji,jj) = qns(ji,jj) - zerp * rcp * sst_m(ji,jj) 
     187                     erp(ji,jj) = zerp 
    190188                  END DO 
    191189               END DO 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90

    r1858 r1859  
    11MODULE trasbc 
    2    !!============================================================================== 
     2   !!====================================================================== 
    33   !!                       ***  MODULE  trasbc  *** 
    44   !! Ocean active tracers:  surface boundary condition 
    5    !!============================================================================== 
    6    !! History :  8.2  !  98-10  (G. Madec, G. Roullet, M. Imbard)  Original code 
    7    !!            8.2  !  01-02  (D. Ludicone)  sea ice and free surface 
    8    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
     5   !!====================================================================== 
     6   !! History :  OPA  !  1998-10  (G. Madec, G. Roullet, M. Imbard)  Original code 
     7   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface 
     8   !!  NEMO      1.0  !  2002-06  (G. Madec)  F90: Free form and module 
     9   !!            3.3  !  2010-05  (Y. Aksenov G. Madec) salt flux + heat associated with emp 
    910   !!---------------------------------------------------------------------- 
    1011 
     
    3132#  include "vectopt_loop_substitute.h90" 
    3233   !!---------------------------------------------------------------------- 
    33    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     34   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    3435   !! $Id$ 
    3536   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    5253      !!         at the surface by evaporation, precipitations and runoff (E-P-R);  
    5354      !!      (3) Fwe, tracer carried with the water that is exchanged.  
     55      !!            - salinity    : salt flux only due to freezing/melting 
     56      !!            sa = sa +  fsalt / rau0 / e3t  for k=1 
    5457      !! 
    5558      !!      Fext, flux through the air-sea interface for temperature and salt:  
     
    7679      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST 
    7780      !!            - salinity    : evaporation, precipitation and runoff 
    78       !!         water has a zero salinity (Fwe=0), thus only Fwi remains: 
    79       !!            sa = sa + emp * sn / e3t   for k=1 
     81      !!         water has a zero salinity  but there is a salt flux due to  
     82      !!         freezing/melting, thus: 
     83      !!            sa = sa + emp * sn / rau0 / e3t   for k=1 
     84      !!                    + fsalt    / rau0 / e3t 
    8085      !!         where emp, the surface freshwater budget (evaporation minus 
    8186      !!         precipitation minus runoff) given in kg/m2/s is divided 
    82       !!         by 1035 kg/m3 (density of ocena water) to obtain m/s.     
     87      !!         by rau0 = 1020 kg/m3 (density of sea water) to obtain m/s.     
    8388      !!         Note: even though Fwe does not appear explicitly for  
    8489      !!         temperature in this routine, the heat carried by the water 
     
    9297      !!         deal with it in this routine. 
    9398      !!           - temperature: Fwe=SST (P-E+R) is added to Fext. 
    94       !!           - salinity:  Fwe = 0, there is no surface flux of salt. 
     99      !!           - salinity:  Fwe sa = sa +  fsalt / rau0 / e3t. 
    95100      !! 
    96101      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated 
     
    103108      INTEGER, INTENT(in) ::   kt     ! ocean time-step index 
    104109      !! 
    105       INTEGER  ::   ji, jj                   ! dummy loop indices 
    106       REAL(wp) ::   zta, zsa, zsrau, zse3t   ! temporary scalars 
     110      INTEGER  ::   ji, jj        ! dummy loop indices 
     111      REAL(wp) ::   z1_e3t_rau0   ! local scalars 
    107112      !!---------------------------------------------------------------------- 
    108113 
     
    113118      ENDIF 
    114119 
    115       zsrau = 1. / rau0             ! initialization 
    116 #if defined key_zco 
    117       zse3t = 1. / e3t_0(1) 
    118 #endif 
    119  
    120120      IF( l_trdtra ) THEN           ! Save ta and sa trends 
    121121         ztrdt(:,:,:) = ta(:,:,:)  
     
    123123      ENDIF 
    124124 
     125!!gm useless staff ??? 
    125126      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
     127!!gm 
    126128 
    127       ! Concentration dillution effect on (t,s) 
    128       DO jj = 2, jpj 
    129          DO ji = fs_2, fs_jpim1   ! vector opt. 
    130 #if ! defined key_zco 
    131             zse3t = 1. / fse3t(ji,jj,1) 
    132 #endif 
    133             IF( lk_vvl) THEN 
    134                zta = r1_rau0_rcp * qns(ji,jj) * zse3t &              ! temperature : heat flux 
    135                 &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux 
    136                zsa = 0.e0                                            ! No salinity concent./dilut. effect 
    137             ELSE 
    138                zta = r1_rau0_rcp * qns(ji,jj) * zse3t                ! temperature : heat flux 
    139                zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect 
    140             ENDIF 
    141             ta(ji,jj,1) = ta(ji,jj,1) + zta                          ! add the trend to the general tracer trend 
    142             sa(ji,jj,1) = sa(ji,jj,1) + zsa 
     129       
     130      IF( lk_vvl ) THEN          ! Variable Volume Layers case   ===>> heat content of mass flux in qns 
     131         DO jj = 2, jpj 
     132            DO ji = fs_2, fs_jpim1   ! vector opt. 
     133               z1_e3t_rau0 = 1./ ( fse3t(ji,jj,1) * rau0 ) 
     134               ta(ji,jj,1) = ta(ji,jj,1) + z1_e3t_rau0 * qns (ji,jj)  * r1_rcp             ! non solar heat flux 
     135               sa(ji,jj,1) = sa(ji,jj,1) + z1_e3t_rau0 * emps(ji,jj)                       ! salt flux (freezing/melting) 
     136            END DO 
    143137         END DO 
    144       END DO 
     138         ! 
     139      ELSE                       ! Constant Volume layers case   ===>> Concentration dillution effect 
     140         DO jj = 2, jpj 
     141            DO ji = fs_2, fs_jpim1   ! vector opt. 
     142               z1_e3t_rau0 = 1./ ( fse3t(ji,jj,1) * rau0 ) 
     143               ta(ji,jj,1) = ta(ji,jj,1) + z1_e3t_rau0 * (  qns (ji,jj) * r1_rcp       &   ! non solar heat flux 
     144                  &                                       + emp (ji,jj) * tn(ji,jj,1)  )   ! concent./dilut. effect 
     145               sa(ji,jj,1) = sa(ji,jj,1) + z1_e3t_rau0 * (  emps(ji,jj)                &   ! salt flux (freezing/melting) 
     146                  &                                       + emp (ji,jj) * sn(ji,jj,1)  )   ! concent./dilut. effect 
     147            END DO 
     148         END DO 
     149      ENDIF 
    145150 
    146151      IF( l_trdtra ) THEN      ! save the sbc trends for diagnostic 
Note: See TracChangeset for help on using the changeset viewer.