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 7351 for branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 – NEMO

Ignore:
Timestamp:
2016-11-28T17:04:10+01:00 (7 years ago)
Author:
emanuelaclementi
Message:

ticket #1805 step 3: /2016/dev_INGV_UKMO_2016 aligned to the trunk at revision 7161

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r5643 r7351  
    1313 
    1414   !!---------------------------------------------------------------------- 
    15    !!   tra_sbc      : update the tracer trend at ocean surface 
    16    !!---------------------------------------------------------------------- 
    17    USE oce             ! ocean dynamics and active tracers 
    18    USE sbc_oce         ! surface boundary condition: ocean 
    19    USE dom_oce         ! ocean space domain variables 
    20    USE phycst          ! physical constant 
    21    USE sbcmod          ! ln_rnf   
    22    USE sbcrnf          ! River runoff   
    23    USE sbcisf          ! Ice shelf    
    24    USE traqsr          ! solar radiation penetration 
    25    USE trd_oce         ! trends: ocean variables 
    26    USE trdtra          ! trends manager: tracers  
     15   !!   tra_sbc       : update the tracer trend at ocean surface 
     16   !!---------------------------------------------------------------------- 
     17   USE oce            ! ocean dynamics and active tracers 
     18   USE sbc_oce        ! surface boundary condition: ocean 
     19   USE dom_oce        ! ocean space domain variables 
     20   USE phycst         ! physical constant 
     21   USE eosbn2         ! Equation Of State 
     22   USE sbcmod         ! ln_rnf   
     23   USE sbcrnf         ! River runoff   
     24   USE sbcisf         ! Ice shelf    
     25   USE iscplini       ! Ice sheet coupling 
     26   USE traqsr         ! solar radiation penetration 
     27   USE trd_oce        ! trends: ocean variables 
     28   USE trdtra         ! trends manager: tracers  
    2729   ! 
    28    USE in_out_manager  ! I/O manager 
    29    USE prtctl          ! Print control 
    30    USE iom 
    31    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    32    USE wrk_nemo        ! Memory Allocation 
    33    USE timing          ! Timing 
    34    USE eosbn2 
     30   USE in_out_manager ! I/O manager 
     31   USE prtctl         ! Print control 
     32   USE iom            ! xIOS server 
     33   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     34   USE wrk_nemo       ! Memory Allocation 
     35   USE timing         ! Timing 
    3536 
    3637   IMPLICIT NONE 
    3738   PRIVATE 
    3839 
    39    PUBLIC   tra_sbc    ! routine called by step.F90 
     40   PUBLIC   tra_sbc   ! routine called by step.F90 
    4041 
    4142   !! * Substitutions 
    42 #  include "domzgr_substitute.h90" 
    4343#  include "vectopt_loop_substitute.h90" 
    4444   !!---------------------------------------------------------------------- 
     
    5757      !!      and add it to the general trend of tracer equations. 
    5858      !! 
    59       !! ** Method : 
    60       !!      Following Roullet and Madec (2000), the air-sea flux can be divided  
    61       !!      into three effects: (1) Fext, external forcing;  
    62       !!      (2) Fwi, concentration/dilution effect due to water exchanged  
    63       !!         at the surface by evaporation, precipitations and runoff (E-P-R);  
    64       !!      (3) Fwe, tracer carried with the water that is exchanged.  
    65       !!            - salinity    : salt flux only due to freezing/melting 
    66       !!            sa = sa +  sfx / rau0 / e3t  for k=1 
     59      !! ** Method :   The (air+ice)-sea flux has two components:  
     60      !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface);  
     61      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice.  
     62      !!               The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe, 
     63      !!             they are simply added to the tracer trend (tsa). 
     64      !!               In linear free surface case (ln_linssh=T), the volume of the 
     65      !!             ocean does not change with the water exchanges at the (air+ice)-sea 
     66      !!             interface. Therefore another term has to be added, to mimic the 
     67      !!             concentration/dilution effect associated with water exchanges. 
    6768      !! 
    68       !!      Fext, flux through the air-sea interface for temperature and salt:  
    69       !!            - temperature : heat flux q (w/m2). If penetrative solar 
    70       !!         radiation q is only the non solar part of the heat flux, the 
    71       !!         solar part is added in traqsr.F routine. 
    72       !!            ta = ta + q /(rau0 rcp e3t)  for k=1 
    73       !!            - salinity    : no salt flux 
    74       !! 
    75       !!      The formulation for Fwb and Fwi vary according to the free  
    76       !!      surface formulation (linear or variable volume).  
    77       !!      * Linear free surface 
    78       !!            The surface freshwater flux modifies the ocean volume 
    79       !!         and thus the concentration of a tracer and the temperature. 
    80       !!         First order of the effect of surface freshwater exchange  
    81       !!         for salinity, it can be neglected on temperature (especially 
    82       !!         as the temperature of precipitations and runoffs is usually 
    83       !!         unknown). 
    84       !!            - temperature : we assume that the temperature of both 
    85       !!         precipitations and runoffs is equal to the SST, thus there 
    86       !!         is no additional flux since in this case, the concentration 
    87       !!         dilution effect is balanced by the net heat flux associated 
    88       !!         to the freshwater exchange (Fwe+Fwi=0): 
    89       !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST 
    90       !!            - salinity    : evaporation, precipitation and runoff 
    91       !!         water has a zero salinity  but there is a salt flux due to  
    92       !!         freezing/melting, thus: 
    93       !!            sa = sa + emp * sn / rau0 / e3t   for k=1 
    94       !!                    + sfx    / rau0 / e3t 
    95       !!         where emp, the surface freshwater budget (evaporation minus 
    96       !!         precipitation minus runoff) given in kg/m2/s is divided 
    97       !!         by rau0 (density of sea water) to obtain m/s.     
    98       !!         Note: even though Fwe does not appear explicitly for  
    99       !!         temperature in this routine, the heat carried by the water 
    100       !!         exchanged through the surface is part of the total heat flux 
    101       !!         forcing and must be taken into account in the global heat 
    102       !!         balance). 
    103       !!      * nonlinear free surface (variable volume, lk_vvl) 
    104       !!         contrary to the linear free surface case, Fwi is properly  
    105       !!         taken into account by using the true layer thicknesses to        
    106       !!         calculate tracer content and advection. There is no need to  
    107       !!         deal with it in this routine. 
    108       !!           - temperature: Fwe=SST (P-E+R) is added to Fext. 
    109       !!           - salinity:  Fwe = 0, there is no surface flux of salt. 
    110       !! 
    111       !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated 
    112       !!                with the tracer surface boundary condition  
    113       !!              - send trends to trdtra module (l_trdtra=T) 
     69      !! ** Action  : - Update tsa with the surface boundary condition trend  
     70      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T) 
    11471      !!---------------------------------------------------------------------- 
    11572      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    116       !! 
    117       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices   
    118       INTEGER  ::   ikt, ikb  
    119       INTEGER  ::   nk_isf 
    120       REAL(wp) ::   zfact, z1_e3t, zdep 
    121       REAL(wp) ::   zalpha, zhk 
     73      ! 
     74      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices   
     75      INTEGER  ::   ikt, ikb              ! local integers 
     76      REAL(wp) ::   zfact, z1_e3t, zdep   ! local scalar 
    12277      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    12378      !!---------------------------------------------------------------------- 
     
    13085         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    13186      ENDIF 
    132  
     87      ! 
    13388      IF( l_trdtra ) THEN                    !* Save ta and sa trends 
    13489         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )  
     
    13691         ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
    13792      ENDIF 
    138  
    139 !!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
     93      ! 
     94!!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist) 
    14095      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration 
    14196         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns 
    142          qsr(:,:) = 0.e0                     ! qsr set to zero 
     97         qsr(:,:) = 0._wp                     ! qsr set to zero 
    14398      ENDIF 
    14499 
     
    146101      !        EMP, SFX and QNS effects 
    147102      !---------------------------------------- 
    148       !                                          Set before sbc tracer content fields 
    149       !                                          ************************************ 
    150       IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1 
    151          !                                         ! ----------------------------------- 
    152          IF( ln_rstart .AND.    &                     ! Restart: read in restart file 
     103      !                             !==  Set before sbc tracer content fields  ==! 
     104      IF( kt == nit000 ) THEN             !* 1st time-step 
     105         IF( ln_rstart .AND.    &               ! Restart: read in restart file 
    153106              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 
    154             IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file' 
     107            IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file' 
    155108            zfact = 0.5_wp 
    156109            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend 
    157110            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend 
    158          ELSE                                         ! No restart or restart not found: Euler forward time stepping 
     111         ELSE                                   ! No restart or restart not found: Euler forward time stepping 
    159112            zfact = 1._wp 
     113            sbc_tsc(:,:,:) = 0._wp 
    160114            sbc_tsc_b(:,:,:) = 0._wp 
    161115         ENDIF 
    162       ELSE                                         ! Swap of forcing fields 
    163          !                                         ! ---------------------- 
     116      ELSE                                !* other time-steps: swap of forcing fields 
    164117         zfact = 0.5_wp 
    165118         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:) 
    166119      ENDIF 
    167       !                                          Compute now sbc tracer content fields 
    168       !                                          ************************************* 
    169  
    170                                                    ! Concentration dilution effect on (t,s) due to   
    171                                                    ! evaporation, precipitation and qns, but not river runoff  
    172                                                 
    173       IF( lk_vvl ) THEN                            ! Variable Volume case  ==>> heat content of mass flux is in qns 
    174          DO jj = 1, jpj 
    175             DO ji = 1, jpi  
    176                sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                              ! non solar heat flux 
    177                sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)                              ! salt flux due to freezing/melting 
     120      !                             !==  Now sbc tracer content fields  ==! 
     121      DO jj = 2, jpj 
     122         DO ji = fs_2, fs_jpim1   ! vector opt. 
     123            sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux 
     124            sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting 
     125         END DO 
     126      END DO 
     127      IF( ln_linssh ) THEN                !* linear free surface   
     128         DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell 
     129            DO ji = fs_2, fs_jpim1   ! vector opt. 
     130               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 
     131               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal) 
    178132            END DO 
    179          END DO 
    180       ELSE                                         ! Constant Volume case ==>> Concentration dilution effect 
     133         END DO                                 !==>> output c./d. term 
     134         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) ) 
     135         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) ) 
     136      ENDIF 
     137      ! 
     138      DO jn = 1, jpts               !==  update tracer trend  ==! 
    181139         DO jj = 2, jpj 
    182             DO ji = fs_2, fs_jpim1   ! vector opt. 
    183                ! temperature : heat flux 
    184                sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                          &   ! non solar heat flux 
    185                   &                  + r1_rau0     * emp(ji,jj)  * tsn(ji,jj,1,jp_tem)       ! concent./dilut. effect 
    186                ! salinity    : salt flux + concent./dilut. effect (both in sfx) 
    187                sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  sfx(ji,jj)                          &   ! salt flux (freezing/melting) 
    188                   &                                + emp(ji,jj) * tsn(ji,jj,1,jp_sal) )      ! concent./dilut. effect 
     140            DO ji = fs_2, fs_jpim1   ! vector opt.   
     141               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t_n(ji,jj,1) 
    189142            END DO 
    190143         END DO 
    191          IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )   ! c/d term on sst 
    192          IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )   ! c/d term on sss 
    193       ENDIF 
    194       ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff   
    195       DO jn = 1, jpts 
    196          DO jj = 2, jpj 
    197             DO ji = fs_2, fs_jpim1   ! vector opt. 
    198                z1_e3t = zfact / fse3t(ji,jj,1) 
    199                tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t 
    200             END DO 
    201          END DO 
    202144      END DO 
    203       !                                          Write in the ocean restart file 
    204       !                                          ******************************* 
    205       IF( lrst_oce ) THEN 
    206          IF(lwp) WRITE(numout,*) 
    207          IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   & 
    208             &                    'at it= ', kt,' date= ', ndastp 
    209          IF(lwp) WRITE(numout,*) '~~~~' 
     145      !                   
     146      IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==! 
    210147         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) ) 
    211148         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) ) 
    212149      ENDIF 
    213150      ! 
    214       ! 
    215151      !---------------------------------------- 
    216152      !       Ice Shelf effects (ISF) 
     
    218154      !---------------------------------------- 
    219155      ! 
    220       IF( nn_isf > 0 ) THEN 
    221          zfact = 0.5e0 
     156!!gm BUG ?   Why no differences between non-linear and linear free surface ? 
     157!!gm         probably taken into account in r1_hisf_tbl : to be verified 
     158      IF( ln_isf ) THEN 
     159         zfact = 0.5_wp 
    222160         DO jj = 2, jpj 
    223161            DO ji = fs_2, fs_jpim1 
    224           
     162               ! 
    225163               ikt = misfkt(ji,jj) 
    226164               ikb = misfkb(ji,jj) 
    227     
     165               ! 
    228166               ! level fully include in the ice shelf boundary layer 
    229                ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst) 
    230167               ! sign - because fwf sign of evapo (rnf sign of precip) 
    231168               DO jk = ikt, ikb - 1 
    232                ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    233169               ! compute trend 
    234                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          & 
    235                      &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) 
    236                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          & 
    237                      &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) 
     170                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                & 
     171                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
     172                     &           * r1_hisf_tbl(ji,jj) 
    238173               END DO 
    239174    
    240175               ! level partially include in ice shelf boundary layer  
    241                ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    242176               ! compute trend 
    243                tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           & 
    244                   &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
    245                tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           & 
    246                   &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)  
     177               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 & 
     178                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
     179                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
     180 
    247181            END DO 
    248182         END DO 
    249183         IF( lrst_oce ) THEN 
    250             IF(lwp) WRITE(numout,*) 
    251             IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   & 
    252                &                    'at it= ', kt,' date= ', ndastp 
    253             IF(lwp) WRITE(numout,*) '~~~~' 
    254             CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)          ) 
     184            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf  (:,:)        ) 
    255185            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) ) 
    256186            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) ) 
     
    269199                  zdep = zfact / h_rnf(ji,jj) 
    270200                  DO jk = 1, nk_rnf(ji,jj) 
    271                                         tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    272                                           &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 
    273                      IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
    274                                           &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep  
     201                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                 & 
     202                                           &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 
     203                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                 & 
     204                                           &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep  
    275205                  END DO 
    276206               ENDIF 
     
    278208         END DO   
    279209      ENDIF 
    280   
    281       IF( l_trdtra )   THEN                      ! send trends for further diagnostics 
     210 
     211      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst 
     212      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss 
     213 
     214      ! 
     215      !---------------------------------------- 
     216      !        Ice Sheet coupling imbalance correction to have conservation 
     217      !---------------------------------------- 
     218      ! 
     219      IF( ln_iscpl .AND. ln_hsb) THEN         ! input of heat and salt due to river runoff  
     220         DO jk = 1,jpk 
     221            DO jj = 2, jpj  
     222               DO ji = fs_2, fs_jpim1 
     223                  zdep = 1._wp / e3t_n(ji,jj,jk)  
     224                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem)                       & 
     225                      &                                         * zdep 
     226                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal)                       & 
     227                      &                                         * zdep   
     228               END DO   
     229            END DO   
     230         END DO 
     231      ENDIF 
     232 
     233      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
    282234         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    283235         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
Note: See TracChangeset for help on using the changeset viewer.