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 2148 for branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA – NEMO

Ignore:
Timestamp:
2010-10-04T15:53:42+02:00 (14 years ago)
Author:
cetlod
Message:

merge LOCEAN 2010 developments branches

Location:
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trabbc.F90

    r2104 r2148  
    3636   REAL(wp) ::   rn_geoflx_cst = 86.4e-3      ! Constant value of geothermal heat flux 
    3737 
    38    INTEGER , DIMENSION(jpi,jpj) ::   nbotlevt   ! ocean bottom level index at T-pt 
    39    REAL(wp), DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend 
     38   INTEGER , DIMENSION(jpi,jpj)         ::   nbotlevt   ! ocean bottom level index at T-pt 
     39   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend 
    4040  
    4141   !! * Substitutions 
  • branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trabbl.F90

    r2123 r2148  
    138138 
    139139 
    140    SUBROUTINE tra_bbl_dif( ptrab, ptraa, kjpt ) 
     140   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt ) 
    141141      !!---------------------------------------------------------------------- 
    142142      !!                  ***  ROUTINE tra_bbl_dif  *** 
     
    155155      !!      convection is satified) 
    156156      !! 
    157       !! ** Action  :   ptraa   increased by the bbl diffusive trend 
     157      !! ** Action  :   pta   increased by the bbl diffusive trend 
    158158      !! 
    159159      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 
     
    161161      !!----------------------------------------------------------------------   
    162162      INTEGER                              , INTENT(in   ) ::   kjpt    ! number of tracers 
    163       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptrab   ! before and now tracer fields 
    164       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptraa   ! tracer trend  
     163      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb   ! before and now tracer fields 
     164      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta   ! tracer trend  
    165165      !! 
    166166      INTEGER  ::   ji, jj, jn   ! dummy loop indices 
    167167      INTEGER  ::   ik           ! local integers 
    168168      REAL(wp) ::   zbtr         ! local scalars 
     169      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! tracer trend  
    169170      !!---------------------------------------------------------------------- 
    170171      ! 
    171       !                                                   ! =========== 
    172172      DO jn = 1, kjpt                                     ! tracer loop 
    173173         !                                                ! =========== 
     174#  if defined key_vectopt_loop 
     175         DO jj = 1, 1   ! vector opt. (forced unrolling) 
     176            DO ji = 1, jpij 
     177#else 
     178         DO jj = 1, jpj 
     179            DO ji = 1, jpi  
     180#endif 
     181               ik = mbkt(ji,jj)                        ! bottom T-level index 
     182               zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S 
     183            END DO 
     184         END DO 
     185         !                                                ! Compute the trend 
    174186#  if defined key_vectopt_loop 
    175187         DO jj = 1, 1   ! vector opt. (forced unrolling) 
     
    181193               ik = mbkt(ji,jj)                            ! bottom T-level index 
    182194               zbtr = e1e2t_r(ji,jj)  / fse3t(ji,jj,ik) 
    183                ptraa(ji,jj,ik,jn) = ptraa(ji,jj,ik,jn)                                                         & 
    184                   &               + (   ahu_bbl(ji,jj) * ( ptrab(ji+1,jj  ,ik,jn) - ptrab(ji  ,jj  ,ik,jn) )   & 
    185                   &                   - ahu_bbl(ji,jj) * ( ptrab(ji  ,jj  ,ik,jn) - ptrab(ji-1,jj  ,ik,jn) )   & 
    186                   &                   + ahv_bbl(ji,jj) * ( ptrab(ji  ,jj+1,ik,jn) - ptrab(ji  ,jj  ,ik,jn) )   & 
    187                   &                   - ahv_bbl(ji,jj) * ( ptrab(ji  ,jj  ,ik,jn) - ptrab(ji  ,jj-1,ik,jn) )   ) * zbtr 
     195               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         & 
     196                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   & 
     197                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   & 
     198                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   & 
     199                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr 
    188200            END DO 
    189201         END DO 
     
    194206    
    195207 
    196    SUBROUTINE tra_bbl_adv( ptrab, ptraa, kjpt ) 
     208   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt ) 
    197209      !!---------------------------------------------------------------------- 
    198210      !!                  ***  ROUTINE trc_bbl  *** 
     
    212224      !!----------------------------------------------------------------------   
    213225      INTEGER                              , INTENT(in   ) ::   kjpt    ! number of tracers 
    214       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptrab   ! before and now tracer fields 
    215       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptraa   ! tracer trend  
     226      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb   ! before and now tracer fields 
     227      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta   ! tracer trend  
    216228      !! 
    217229      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
     
    240252                  !                                               ! up  -slope T-point (shelf bottom point) 
    241253                  zbtr = e1e2t_r(iis,jj) / fse3t(iis,jj,ikus) 
    242                   ztra = zu_bbl * ( ptrab(iid,jj,ikus,jn) - ptrab(iis,jj,ikus,jn) ) * zbtr 
    243                   ptraa(iis,jj,ikus,jn) = ptraa(iis,jj,ikus,jn) + ztra 
     254                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 
     255                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 
    244256                  !                    
    245257                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column) 
    246258                     zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,jk) 
    247                      ztra = zu_bbl * ( ptrab(iid,jj,jk+1,jn) - ptrab(iid,jj,jk,jn) ) * zbtr 
    248                      ptraa(iid,jj,jk,jn) = ptraa(iid,jj,jk,jn) + ztra 
     259                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 
     260                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 
    249261                  END DO 
    250262                  !  
    251263                  zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,ikud) 
    252                   ztra = zu_bbl * ( ptrab(iis,jj,ikus,jn) - ptrab(iid,jj,ikud,jn) ) * zbtr 
    253                   ptraa(iid,jj,ikud,jn) = ptraa(iid,jj,ikud,jn) + ztra 
     264                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 
     265                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra 
    254266               ENDIF 
    255267               ! 
     
    262274                  ! up  -slope T-point (shelf bottom point) 
    263275                  zbtr = e1e2t_r(ji,ijs) / fse3t(ji,ijs,ikvs) 
    264                   ztra = zv_bbl * ( ptrab(ji,ijd,ikvs,jn) - ptrab(ji,ijs,ikvs,jn) ) * zbtr 
    265                   ptraa(ji,ijs,ikvs,jn) = ptraa(ji,ijs,ikvs,jn) + ztra 
     276                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 
     277                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 
    266278                  !                    
    267279                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column) 
    268280                     zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,jk) 
    269                      ztra = zv_bbl * ( ptrab(ji,ijd,jk+1,jn) - ptrab(ji,ijd,jk,jn) ) * zbtr 
    270                      ptraa(ji,ijd,jk,jn) = ptraa(ji,ijd,jk,jn)  + ztra 
     281                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 
     282                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra 
    271283                  END DO 
    272284                  !                                               ! down-slope T-point (deep bottom point) 
    273285                  zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,ikvd) 
    274                   ztra = zv_bbl * ( ptrab(ji,ijs,ikvs,jn) - ptrab(ji,ijd,ikvd,jn) ) * zbtr 
    275                   ptraa(ji,ijd,ikvd,jn) = ptraa(ji,ijd,ikvd,jn) + ztra 
     286                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 
     287                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra 
    276288               ENDIF 
    277289            END DO 
     
    370382#endif 
    371383            ik = mbkt(ji,jj)                        ! bottom T-level index 
    372             ztb (ji,jj) = tsb(ji,jj,ik,jp_tem)      ! bottom before T and S 
    373             zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) 
     384            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S 
     385            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1) 
    374386            zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth 
    375387            ! 
     
    575587      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1) 
    576588 
     589 
    577590      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl 
    578591         ! 
  • branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/tranxt.F90

    r2120 r2148  
    1515   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf 
    1616   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option 
    17    !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
     17   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA 
     18   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
    1819   !!---------------------------------------------------------------------- 
    1920 
    2021   !!---------------------------------------------------------------------- 
    21    !!   tra_nxt       : time stepping on temperature and salinity 
    22    !!   tra_nxt_fix   : time stepping on temperature and salinity : fixed    volume case 
    23    !!   tra_nxt_vvl   : time stepping on temperature and salinity : variable volume case 
     22   !!   tra_nxt       : time stepping on tracers 
     23   !!   tra_nxt_fix   : time stepping on tracers : fixed    volume case 
     24   !!   tra_nxt_vvl   : time stepping on tracers : variable volume case 
    2425   !!---------------------------------------------------------------------- 
    2526   USE oce             ! ocean dynamics and tracers variables 
    2627   USE dom_oce         ! ocean space and time domain variables  
     28   USE sbc_oce         ! surface boundary condition: ocean 
    2729   USE zdf_oce         ! ??? 
    2830   USE domvvl          ! variable volume 
     
    3840   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3941   USE prtctl          ! Print control 
     42   USE traqsr          ! penetrative solar radiation (needed for nksr) 
    4043   USE traswp          ! swap array 
    4144   USE agrif_opa_update 
     
    4952   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt 
    5053 
     54   REAL(wp)                 ::   rbcp            ! Brown & Campana parameters for semi-implicit hpg 
    5155   REAL(wp), DIMENSION(jpk) ::   r2dt   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 
    5256 
     
    96100         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 
    97101         IF(lwp) WRITE(numout,*) '~~~~~~~' 
     102         ! 
     103         rbcp    = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)       ! Brown & Campana parameter for semi-implicit hpg 
    98104      ENDIF 
    99105 
     
    107113#endif 
    108114 
    109 #if defined key_obc 
     115#if defined key_obc  
    110116      IF( lk_obc )   CALL obc_tra( kt )  ! OBC open boundaries 
    111117#endif 
     
    114120#endif 
    115121#if defined key_agrif 
    116       CALL Agrif_tra                     ! AGRIF zoom boundaries 
     122      CALL Agrif_tra                   ! AGRIF zoom boundaries 
    117123#endif 
    118124 
     
    133139 
    134140      ! Leap-Frog + Asselin filter time stepping 
    135       IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, tsb, tsn, tsa, jpts )  ! variable volume level (vvl)      
    136       ELSE                  ;   CALL tra_nxt_fix( kt, tsb, tsn, tsa, jpts )  ! fixed    volume level  
     141      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)      
     142      ELSE                  ;   CALL tra_nxt_fix( kt,        tsb, tsn, tsa, jpts )  ! fixed    volume level  
    137143      ENDIF 
    138144 
     
    176182      !!              - swap tracer fields to prepare the next time_step. 
    177183      !!                This can be summurized for tempearture as: 
    178       !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T 
    179       !!             ztm = 0                   otherwise 
     184      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T 
     185      !!             ztm = 0                                   otherwise 
     186      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp) 
    180187      !!             tb  = tn + atfp*[ tb - 2 tn + ta ] 
    181       !!             tn  = ta  
     188      !!             tn  = ta   
    182189      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call) 
    183190      !! 
     
    185192      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    186193      !!---------------------------------------------------------------------- 
    187       INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index 
    188       INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers 
    189       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields 
    190       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields 
    191       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend 
     194      INTEGER , INTENT(in   )                               ::  kt       ! ocean time-step index 
     195      INTEGER , INTENT(in   )                               ::  kjpt     ! number of tracers 
     196      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields 
     197      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields 
     198      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend 
    192199      !! 
    193200      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices 
    194       REAL(wp) :: ztm, ztf     ! temporary scalars 
     201      REAL(wp) :: ztd, ztm         ! temporary scalars 
    195202      !!---------------------------------------------------------------------- 
    196203 
     
    205212         !                                             ! (only swap) 
    206213         DO jn = 1, kjpt 
    207             DO jk = 1, jpkm1                               
     214            DO jk = 1, jpkm1 
    208215               ptn(:,:,jk,jn) = pta(:,:,jk,jn)     ! ptb <-- ptn 
    209216            END DO 
     
    219226                  DO jj = 1, jpj 
    220227                     DO ji = 1, jpi 
    221                         ztm = 0.25 * ( pta(ji,jj,jk,jn) + 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )  ! mean pt 
    222                         ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )  ! Asselin filter on pt  
    223                         ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                      ! ptb <-- filtered ptn  
    224                         ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                            ! ptn <-- pta 
    225                         pta(ji,jj,jk,jn) = ztm                                                           ! pta <-- mean pt 
     228                        ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    !  time laplacian on tracers 
     229                        ztm = ptn(ji,jj,jk,jn) + rbcp * ztd                                 !  used for both Asselin and Brown & Campana filters 
     230                        ! 
     231                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn  
     232                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta 
     233                        pta(ji,jj,jk,jn) = ztm                                              ! pta <-- Brown & Campana average 
    226234                     END DO 
    227235                  END DO 
     
    235243                  DO jj = 1, jpj 
    236244                     DO ji = 1, jpi 
    237                         ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )       ! Asselin filter on t  
    238                         ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                     ! ptb <-- filtered ptn  
    239                         ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                           ! ptn <-- pta 
     245                        ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    ! time laplacian on tracers 
     246                        !  
     247                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn  
     248                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta 
    240249                     END DO 
    241250                  END DO 
    242251               END DO 
    243252            END DO 
    244             ! 
    245253         ENDIF 
    246254         ! 
     
    250258 
    251259 
    252    SUBROUTINE tra_nxt_vvl( kt, ptb, ptn, pta, kjpt ) 
     260   SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt ) 
    253261      !!---------------------------------------------------------------------- 
    254262      !!                   ***  ROUTINE tra_nxt_vvl  *** 
     
    263271      !!              - swap tracer fields to prepare the next time_step. 
    264272      !!                This can be summurized for tempearture as: 
    265       !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T 
    266       !!                  /(e3t_a   +2*e3t_n   +e3t_b   )     
    267       !!             ztm = 0                                otherwise 
     273      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T 
     274      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )    
     275      !!             ztm = 0                                                       otherwise 
    268276      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
    269277      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] ) 
     
    274282      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    275283      !!---------------------------------------------------------------------- 
    276       INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index 
    277       INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers 
    278       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields 
    279       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields 
    280       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend 
     284      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index 
     285      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator) 
     286      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers 
     287      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields 
     288      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields 
     289      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend 
    281290      !!      
    282       INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices 
    283       REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar 
    284       REAL(wp) ::   ze3mr, ze3fr                           !    -         - 
    285       REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         - 
     291      INTEGER  ::   ji, jj, jk, jn                 ! dummy loop indices 
     292      REAL(wp) ::   ztc_a , ztc_n , ztc_b          ! temporary scalar 
     293      REAL(wp) ::   ztc_f , ztc_d , ztc_m          !    -         - 
     294      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a         !    -         - 
     295      REAL(wp) ::   ze3t_f, ze3t_d, ze3t_m         !    -         - 
     296      REAL     ::   zfact1, zfact2                 !    -         - 
    286297      !!---------------------------------------------------------------------- 
    287298 
     
    293304      ! 
    294305      ! 
    295       IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
    296          !                                             ! (only swap) 
     306      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step 
     307         !                                           ! (only swap) 
    297308         DO jn = 1, kjpt 
    298             DO jk = 1, jpkm1                               
    299                ptn(:,:,jk,jn) = pta(:,:,jk,jn)         ! ptb <-- ptn 
     309            DO jk = 1, jpkm1 
     310               ptn(:,:,jk,jn) = pta(:,:,jk,jn)       ! ptb <-- ptn 
    300311            END DO 
    301312         END DO 
    302313         !                                               
    303       ELSE                                              ! general case (Leapfrog + Asselin filter 
     314      ELSE                                           ! general case (Leapfrog + Asselin filter) 
    304315         ! 
    305          !                                              ! ----------------------- ! 
    306          IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case ! 
    307             !                                           ! ----------------------- ! 
    308             DO jn = 1, kjpt                                
     316         !                                           ! ----------------------- ! 
     317         IF( ln_dynhpg_imp ) THEN                    ! semi-implicite hpg case ! 
     318            !                                        ! ----------------------- ! 
     319            DO jn = 1, kjpt                           
    309320               DO jk = 1, jpkm1 
    310321                  DO jj = 1, jpj 
     
    314325                        ze3t_a = fse3t_a(ji,jj,jk) 
    315326                        !                                         ! tracer content at Before, now and after 
    316                         ztcb = ptb(ji,jj,jk,jn) *  ze3t_b   
    317                         ztcn = ptn(ji,jj,jk,jn) *  ze3t_n   
    318                         ztca = pta(ji,jj,jk,jn) *  ze3t_a   
    319                         ! 
    320                         !                                         ! Asselin filter on thickness and tracer content 
    321                         ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 
    322                         ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   )  
    323                         ! 
    324                         !                                         ! filtered tracer including the correction  
    325                         ze3fr = 1.e0  / ( ze3t_n + ze3t_f ) 
    326                         ztf   = ze3fr * ( ztcn   + ztc_f  ) 
    327                         !                                         ! mean thickness and tracer 
    328                         ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b ) 
    329                         ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   ) 
    330                         !!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !! 
    331                         !!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr 
     327                        ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b   
     328                        ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n   
     329                        ztc_a  = pta(ji,jj,jk,jn) * ze3t_a   
     330                        !                                         ! Time laplacian on tracer contents 
     331                        !                                         ! used for both Asselin and Brown & Campana filters 
     332                        ze3t_d =  ze3t_b - 2.* ze3t_n + ze3t_a  
     333                        ztc_d  =  ztc_b  - 2.* ztc_n  + ztc_a   
     334                        !                                         ! Asselin Filter on thicknesses and tracer contents 
     335                        ztc_f  = ztc_n  + atfp * ztc_d 
     336                        ztc_m  = ztc_n  + rbcp * ztc_d 
     337                        !                                          
     338                        ze3t_f = 1.0 / ( ze3t_n + atfp * ze3t_d ) 
     339                        ze3t_m = 1.0 / ( ze3t_n + rbcp * ze3t_d ) 
    332340                        !                                         ! swap of arrays 
    333                         ptb(ji,jj,jk,jn) = ztf                             ! ptb <-- ptn + filter 
    334                         ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)              ! ptn <-- pta 
    335                         pta(ji,jj,jk,jn) = ztm                             ! pta <-- mean t 
     341                        ptb(ji,jj,jk,jn) = ztc_f * ze3t_f              ! ptb <-- ptn filtered 
     342                        pta(ji,jj,jk,jn) = ztc_m * ze3t_m              ! pta <-- Brown & Campana average 
     343                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)            ! ptn <-- pta 
    336344                     END DO 
    337345                  END DO 
    338346               END DO 
    339347            END DO 
    340             !                                        ! ----------------------- ! 
    341          ELSE                                        !    explicit hpg case    ! 
    342             !                                        ! ----------------------- ! 
    343             DO jn = 1, kjpt       
    344                DO jk = 1, jpkm1 
    345                   DO jj = 1, jpj 
    346                      DO ji = 1, jpi 
    347                         ze3t_b = fse3t_b(ji,jj,jk) 
    348                         ze3t_n = fse3t_n(ji,jj,jk) 
    349                         ze3t_a = fse3t_a(ji,jj,jk) 
    350                         !                                         ! tracer content at Before, now and after 
    351                         ztcb = ptb(ji,jj,jk,jn) *  ze3t_b   
    352                         ztcn = ptn(ji,jj,jk,jn) *  ze3t_n    
    353                         ztca = pta(ji,jj,jk,jn) *  ze3t_a   
    354                         ! 
    355                         !                                         ! Asselin filter on thickness and tracer content 
    356                         ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 
    357                         ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   )  
    358                         ! 
    359                         !                                         ! filtered tracer including the correction  
    360                         ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 
    361                         ztf   =  ( ztcn  + ztc_f ) * ze3fr 
    362                         !                                         ! swap of arrays 
    363                         ptb(ji,jj,jk,jn) = ztf                    ! tb <-- tn filtered 
    364                         ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)       ! tn <-- ta 
     348            !                                           ! ----------------------- ! 
     349         ELSE                                           !    explicit hpg case    ! 
     350            !                                           ! ----------------------- ! 
     351            IF( cdtype == 'TRA' ) THEN 
     352               !               
     353               DO jn = 1, kjpt       
     354                  DO jk = 1, jpkm1 
     355                     zfact1 = atfp * rdttra(jk) 
     356                     zfact2 = zfact1 / rau0 
     357                     DO jj = 1, jpj 
     358                        DO ji = 1, jpi 
     359                           ze3t_b = fse3t_b(ji,jj,jk) 
     360                           ze3t_n = fse3t_n(ji,jj,jk) 
     361                           ze3t_a = fse3t_a(ji,jj,jk) 
     362                           !                                         ! tracer content at Before, now and after 
     363                           ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
     364                           ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n 
     365                           ztc_a  = pta(ji,jj,jk,jn) * ze3t_a 
     366                           ! 
     367                           ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a ) 
     368                           ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b ) 
     369 
     370                           IF( jk == 1 ) THEN 
     371                               ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 
     372                               IF( jn == jp_tem )   ztc_f  = ztc_f  - zfact1 * ( sbc_hc_n(ji,jj) - sbc_hc_b(ji,jj) ) 
     373                               IF( jn == jp_sal )   ztc_f  = ztc_f  - zfact1 * ( sbc_sc_n(ji,jj) - sbc_sc_b(ji,jj) ) 
     374                           ENDIF 
     375                           IF( jn == jp_tem .AND. ln_traqsr .AND. jk <= nksr )  & 
     376                           &                        ztc_f  = ztc_f  - zfact1 * ( qsr_hc_n(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )  
     377 
     378                           ze3t_f = 1.e0 / ze3t_f 
     379                           ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered 
     380                           ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta 
     381                        END DO 
    365382                     END DO 
    366383                  END DO 
    367384               END DO 
    368             END DO 
     385               ! 
     386            ELSE IF( cdtype == 'TRC' ) THEN 
     387               !               
     388               DO jn = 1, kjpt       
     389                  DO jk = 1, jpkm1 
     390                     DO jj = 1, jpj 
     391                        DO ji = 1, jpi 
     392                           ze3t_b = fse3t_b(ji,jj,jk) 
     393                           ze3t_n = fse3t_n(ji,jj,jk) 
     394                           ze3t_a = fse3t_a(ji,jj,jk) 
     395                           !                                         ! tracer content at Before, now and after 
     396                           ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
     397                           ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n 
     398                           ztc_a  = pta(ji,jj,jk,jn) * ze3t_a 
     399                           ! 
     400                           ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a ) 
     401                           ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b  ) 
     402 
     403                           ze3t_f = 1.e0 / ze3t_f 
     404                           ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered 
     405                           ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta 
     406                        END DO 
     407                     END DO 
     408                  END DO 
     409               END DO 
     410               ! 
     411            END IF 
    369412            ! 
    370413         ENDIF 
  • branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/traqsr.F90

    r2024 r2148  
    2727   USE iom             ! I/O manager 
    2828   USE fldread         ! read input fields 
     29   USE restart         ! ocean restart 
    2930 
    3031   IMPLICIT NONE 
     
    4748   ! Module variables 
    4849   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
    49    INTEGER ::   nksr   ! levels below which the light cannot penetrate ( depth larger than 391 m) 
     50   INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m) 
    5051   REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption 
    5152 
     
    9596      REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars 
    9697      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         - 
     98      REAL(wp) ::   z1_e3t, zfact        !    -         - 
    9799      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace 
    98100      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace 
     
    111113         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = 0. 
    112114      ENDIF 
     115 
     116      !                                        Set before qsr tracer content field 
     117      !                                        *********************************** 
     118      IF( kt == nit000 ) THEN                     ! Set the forcing field at nit000 - 1 
     119         !                                        ! ----------------------------------- 
     120         IF( ln_rstart .AND.    &                    ! Restart: read in restart file 
     121              & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN 
     122            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file' 
     123            zfact = 0.5e0 
     124            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux 
     125         ELSE                                           ! No restart or restart not found: Euler forward time stepping 
     126            zfact = 1.e0 
     127            qsr_hc_b(:,:,:) = 0.e0 
     128         ENDIF 
     129      ELSE                                        ! Swap of forcing field 
     130         !                                        ! --------------------- 
     131         zfact = 0.5e0 
     132         qsr_hc_b(:,:,:) = qsr_hc_n(:,:,:) 
     133      ENDIF 
     134      !                                        Compute now qsr tracer content field 
     135      !                                        ************************************ 
    113136       
    114137      !                                           ! ============================================== ! 
     
    118141            DO jj = 2, jpjm1 
    119142               DO ji = fs_2, fs_jpim1   ! vector opt. 
    120                   ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)  
     143                  qsr_hc_n(ji,jj,jk) = ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)  
    121144               END DO 
    122145            END DO 
     
    175198               ! 
    176199               DO jk = 1, nksr                                        ! compute and add qsr trend to ta 
    177                   tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
     200                  qsr_hc_n(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
    178201               END DO 
    179202               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero 
     
    182205            ELSE                                                 !*  Constant Chlorophyll 
    183206               DO jk = 1, nksr 
    184                   tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + etot3(:,:,jk) * qsr(:,:) 
     207                  qsr_hc_n(:,:,jk) = etot3(:,:,jk) * qsr(:,:) 
    185208               END DO 
    186209            ENDIF 
     
    194217               DO jj = 2, jpjm1 
    195218                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    196                      tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + etot3(ji,jj,jk) * qsr(ji,jj) 
     219                     qsr_hc_n(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) 
    197220                  END DO 
    198221               END DO 
     
    200223            ! 
    201224         ENDIF 
     225         ! 
     226      ENDIF 
     227      !                                        Add to the general trend 
     228      !                                        ************************ 
     229      DO jk = 1, nksr 
     230         DO jj = 2, jpjm1  
     231            DO ji = fs_2, fs_jpim1   ! vector opt. 
     232               z1_e3t = zfact / fse3t(ji,jj,jk) 
     233               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc_n(ji,jj,jk) ) * z1_e3t 
     234            END DO 
     235         END DO 
     236      END DO 
     237      ! 
     238      IF( lrst_oce ) THEN   !                  Write in the ocean restart file 
     239         !                                     ******************************* 
     240         IF(lwp) WRITE(numout,*) 
     241         IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ',   & 
     242            &                    'at it= ', kt,' date= ', ndastp 
     243         IF(lwp) WRITE(numout,*) '~~~~' 
     244         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc_n ) 
    202245         ! 
    203246      ENDIF 
  • branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trasbc.F90

    r2113 r2148  
    77   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface 
    88   !!  NEMO      1.0  !  2002-06  (G. Madec)  F90: Free form and module 
    9    !!            3.3  !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
     9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
     10   !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    2223   USE in_out_manager  ! I/O manager 
    2324   USE prtctl          ! Print control 
     25   USE restart         ! ocean restart 
    2426   USE sbcrnf          ! River runoff   
    2527   USE sbcmod          ! ln_rnf   
     28   USE iom 
    2629 
    2730   IMPLICIT NONE 
     
    103106      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    104107      !! 
    105       INTEGER  ::   ji, jj, jk      ! dummy loop indices   
    106       REAL(wp) ::   zta, zsa        ! local scalars, adjustment to temperature and salinity   
    107       REAL(wp) ::   zata, zasa      ! local scalars, calculations of automatic change to temp & sal due to vvl (done elsewhere)   
     108      INTEGER  ::   ji, jj, jk           ! dummy loop indices   
     109      REAL(wp) ::   zta, zsa, zrnf       ! local scalars, adjustment to temperature and salinity   
    108110      REAL(wp) ::   zsrau, zse3t, zdep   ! local scalars, 1/density, 1/height of box, 1/height of effected water column   
    109111      REAL(wp) ::   zdheat, zdsalt       ! total change of temperature and salinity   
     112      REAL(wp) ::   zfact, z1_e3t        !  
    110113      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds 
    111114      !!---------------------------------------------------------------------- 
     
    118121 
    119122      zsrau = 1. / rau0             ! initialization 
    120 #if defined key_zco 
    121       zse3t = 1. / e3t_0(1) 
    122 #endif 
    123123 
    124124      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
     
    127127      ENDIF 
    128128 
    129       IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
    130  
    131       ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff   
     129!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
     130      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration 
     131         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns 
     132         qsr(:,:) = 0.e0                     ! qsr set to zero 
     133      ENDIF 
     134 
     135      !                                          Set before sbc tracer content fields 
     136      !                                          ************************************ 
     137      IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1 
     138         !                                         ! ----------------------------------- 
     139         IF( ln_rstart .AND.    &                     ! Restart: read in restart file 
     140              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 
     141            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file' 
     142            zfact = 0.5e0 
     143            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_hc_b )   ! before heat content sbc trend 
     144            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_sc_b )   ! before salt content sbc trend 
     145         ELSE                                         ! No restart or restart not found: Euler forward time stepping 
     146            zfact = 1.e0 
     147            sbc_hc_b(:,:) = 0.e0 
     148            sbc_sc_b(:,:) = 0.e0 
     149         ENDIF 
     150      ELSE                                         ! Swap of forcing fields 
     151         !                                         ! ---------------------- 
     152         zfact = 0.5e0 
     153         sbc_hc_b(:,:) = sbc_hc_n(:,:) 
     154         sbc_sc_b(:,:) = sbc_sc_n(:,:) 
     155      ENDIF 
     156      !                                          Compute now sbc tracer content fields 
     157      !                                          ************************************* 
     158 
     159                                                   ! Concentration dilution effect on (t,s) due to   
     160                                                   ! evaporation, precipitation and qns, but not river runoff  
     161                                                
     162      IF( lk_vvl ) THEN                            ! Variable Volume case 
     163         DO jj = 2, jpj 
     164            DO ji = fs_2, fs_jpim1   ! vector opt. 
     165               ! temperature : heat flux + cooling/heating effet of EMP flux 
     166               sbc_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 
     167               ! concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration 
     168               sbc_sc_n(ji,jj) = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal) 
     169            END DO 
     170         END DO 
     171      ELSE                                         ! Constant Volume case 
     172         DO jj = 2, jpj 
     173            DO ji = fs_2, fs_jpim1   ! vector opt. 
     174               ! temperature : heat flux 
     175               sbc_hc_n(ji,jj) = ro0cpr * qns(ji,jj) 
     176               ! salinity    : salt flux + concent./dilut. effect (both in emps) 
     177               sbc_sc_n(ji,jj) = zsrau * emps(ji,jj) * tsn(ji,jj,1,jp_sal) 
     178            END DO 
     179         END DO 
     180      ENDIF 
     181                                                   ! Concentration dilution effect on (t,s) due to   
     182                                                   ! river runoff without T, S and depth attributes 
     183      IF( ln_rnf ) THEN  
     184         ! 
     185         IF( lk_vvl ) THEN                            ! Variable Volume case  
     186            !                                         ! cooling/heating effect of runoff & No salinity concent./dilut. effect 
     187            DO jj = 2, jpj 
     188               DO ji = fs_2, fs_jpim1   ! vector opt. 
     189                  sbc_hc_n(ji,jj) = sbc_hc_n(ji,jj) + zsrau * rnf(ji,jj) * tsn(ji,jj,1,jp_tem) 
     190               END DO 
     191            END DO 
     192         ELSE                                         ! Constant Volume case 
     193            !                                         ! concent./dilut. effect only 
     194            DO jj = 2, jpj  
     195               DO ji = fs_2, fs_jpim1   ! vector opt. 
     196                  sbc_sc_n(ji,jj) = sbc_sc_n(ji,jj) - zsrau * rnf(ji,jj) * tsn(ji,jj,1,jp_sal) 
     197               END DO 
     198            END DO 
     199         ENDIF 
     200         !   
     201      ENDIF 
     202      !                                          Add to the general trend 
     203      !                                          ************************ 
    132204      DO jj = 2, jpj 
    133205         DO ji = fs_2, fs_jpim1   ! vector opt. 
    134 #if ! defined key_zco 
    135             zse3t = 1. / fse3t(ji,jj,1) 
    136 #endif 
    137             IF( lk_vvl) THEN 
    138                ! temperature : heat flux and heat content of EMP flux 
    139                zta = ( ro0cpr * qns(ji,jj) - emp(ji,jj) * zsrau * tsn(ji,jj,1,jp_tem) ) * zse3t 
    140                ! Salinity    : concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration 
    141                zsa = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t 
    142             ELSE 
    143                zta =  ro0cpr * qns(ji,jj) * zse3t                         ! temperature : heat flux  
    144                zsa =  emps(ji,jj) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t   ! salinity :  concent./dilut. effect  
    145             ENDIF 
    146             tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta                  ! add the trend to the general tracer trend 
    147             tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa 
     206            z1_e3t = zfact / fse3t(ji,jj,1) 
     207            tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + ( sbc_hc_b(ji,jj) + sbc_hc_n(ji,jj) ) * z1_e3t 
     208            tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + ( sbc_sc_b(ji,jj) + sbc_sc_n(ji,jj) ) * z1_e3t 
    148209         END DO 
    149210      END DO 
     211      !                                          Write in the ocean restart file 
     212      !                                          ******************************* 
     213      IF( lrst_oce ) THEN 
     214         IF(lwp) WRITE(numout,*) 
     215         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   & 
     216            &                    'at it= ', kt,' date= ', ndastp 
     217         IF(lwp) WRITE(numout,*) '~~~~' 
     218         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_hc_n ) 
     219         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_sc_n ) 
     220      ENDIF 
    150221 
    151222      IF( ln_rnf .AND. ln_rnf_att ) THEN        ! Concentration / dilution effect on (t,s) due to river runoff   
     223        ! 
    152224        DO jj = 1, jpj   
    153225           DO ji = 1, jpi   
    154               zdep = 1. / rnf_dep(ji,jj)   
    155               zse3t= 1. / fse3t(ji,jj,1)   
    156               IF( rnf_tem(ji,jj) == -999 )   rnf_tem(ji,jj) = tsn(ji,jj,1,jp_tem)   ! if not specified set runoff temp to be sst   
     226              zdep  = 1. / rnf_dep(ji,jj)   
     227              zse3t = 1. / fse3t(ji,jj,1)   
     228              rnf_dep(ji,jj) = 0.e0 
     229              DO jk = 1, rnf_mod_dep(ji,jj)                          ! recalculates rnf_dep to be the depth   
     230                rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk)    ! in metres to the bottom of the relevant grid box   
     231              END DO 
     232              IF( rnf_tmp(ji,jj) == -999 )   rnf_tmp(ji,jj) = tsn(ji,jj,1,jp_tem)   ! if not specified set runoff temp to be sst   
    157233   
    158234              IF( rnf(ji,jj) > 0.e0 ) THEN   
    159    
     235                ! 
     236                zrnf = rnf(ji,jj) * zsrau * zdep 
    160237                IF( lk_vvl ) THEN   
    161238                  ! indirect flux, concentration or dilution effect : force a dilution effect in all levels  
     
    163240                  zdsalt = 0.e0   
    164241                  DO jk = 1, rnf_mod_dep(ji,jj)   
    165                     zta = -tsn(ji,jj,jk,jp_tem) * rnf(ji,jj) * zsrau * zdep   
    166                     zsa = -tsn(ji,jj,jk,jp_sal) * rnf(ji,jj) * zsrau * zdep   
     242                    zta = -tsn(ji,jj,jk,jp_tem) * zrnf   
     243                    zsa = -tsn(ji,jj,jk,jp_sal) * zrnf   
    167244                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend 
    168245                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa 
     
    171248                  END DO   
    172249                  ! negate this total change in heat and salt content from top level    !!gm I don't understand this 
    173                   zta = -zdheat * zse3t   
    174                   zsa = -zdsalt * zse3t   
    175                   tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta            ! add the trend to the general tracer trend 
    176                   tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa 
    177      
    178                   ! direct flux   
    179                   zta = rnf_tem(ji,jj) * rnf(ji,jj) * zsrau * zdep   
    180                   zsa = rnf_sal(ji,jj) * rnf(ji,jj) * zsrau * zdep   
     250                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) - zdheat * zse3t 
     251                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) - zdsalt * zse3t 
    181252     
    182253                  DO jk = 1, rnf_mod_dep(ji,jj)   
    183                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend 
    184                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa 
     254                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + rnf_tmp(ji,jj) * zrnf 
     255                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + rnf_sal(ji,jj) * zrnf 
    185256                  END DO   
    186257                ELSE   
    187258                  DO jk = 1, rnf_mod_dep(ji,jj)   
    188                     zta = ( rnf_tem(ji,jj) - tsn(ji,jj,jk,jp_tem) ) * rnf(ji,jj) * zsrau * zdep   
    189                     zsa = ( rnf_sal(ji,jj) - tsn(ji,jj,jk,jp_sal) ) * rnf(ji,jj) * zsrau * zdep   
    190                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend 
    191                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa 
     259                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( rnf_tmp(ji,jj) - tsn(ji,jj,jk,jp_tem) ) * zrnf 
     260                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + ( rnf_sal(ji,jj) - tsn(ji,jj,jk,jp_sal) ) * zrnf 
    192261                  END DO   
    193262                ENDIF   
    194263   
    195               ELSEIF( rnf(ji,jj) < 0.e0) THEN   ! for use in baltic when flow is out of domain, want no change in temp and sal   
     264              ELSEIF( rnf(ji,jj) < 0.e0 ) THEN   ! for use in baltic when flow is out of domain, want no change in temp and sal   
    196265   
    197266                IF( lk_vvl ) THEN   
    198267                  ! calculate automatic adjustment to sal and temp due to dilution/concentraion effect    
    199                   zata = tsn(ji,jj,1,jp_tem) * rnf(ji,jj) * zsrau * zse3t   
    200                   zasa = tsn(ji,jj,1,jp_sal) * rnf(ji,jj) * zsrau * zse3t   
    201                   tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zata                  ! add the trend to the general tracer trend 
    202                   tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zasa 
     268                  zrnf = rnf(ji,jj) * zsrau * zse3t 
     269                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * zrnf  
     270                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * zrnf 
    203271                ENDIF   
    204272   
     
    207275           END DO   
    208276        END DO   
    209  
    210       ELSE IF( ln_rnf ) THEN      ! Concentration dilution effect on (t,s) due to runoff without T, S and depth attributes 
    211  
    212  
    213         DO jj = 2, jpj 
    214            DO ji = fs_2, fs_jpim1   ! vector opt. 
    215 #if ! defined key_zco 
    216               zse3t = 1. / fse3t(ji,jj,1) 
    217 #endif 
    218               IF( lk_vvl) THEN 
    219                  zta =    rnf(ji,jj) * zsrau * tsn(ji,jj,1,jp_tem) * zse3t       ! & cooling/heating effect of runoff 
    220                  zsa =    0.e0                                                   ! No salinity concent./dilut. effect 
    221               ELSE 
    222                  zta =    0.0                                            ! temperature : heat flux  
    223                  zsa =  - rnf(ji,jj) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t     ! salinity :  concent./dilut. effect 
    224               ENDIF 
    225               tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta          ! add the trend to the general tracer trend 
    226               tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa 
    227            END DO 
    228         END DO 
    229   
     277        ! 
    230278      ENDIF   
    231279 
Note: See TracChangeset for help on using the changeset viewer.