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 12590 – NEMO

Changeset 12590


Ignore:
Timestamp:
2020-03-23T22:16:19+01:00 (3 years ago)
Author:
techene
Message:

all: add e3 substitute, OCE/DOM/domzgr_substitute.h90: correct a bug for e3f

Location:
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE
Files:
23 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90

    r12583 r12590  
    1515#   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1.+r3u(i,j,t)*umask(i,j,k))) 
    1616#   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1.+r3v(i,j,t)*vmask(i,j,k))) 
    17 #   define  e3f(i,j,k,t)   (e3f_0(i,j,k)*(1.+r3f(i,j,t)*fmask(i,j,k))) 
     17#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1.+r3f(i,j)*fmask(i,j,k))) 
    1818#   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1.+r3t(i,j,t))) 
    1919#   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1.+r3u(i,j,t))) 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/divhor.F90

    r12377 r12590  
    2121   USE dom_oce         ! ocean space and time domain 
    2222   USE sbc_oce, ONLY : ln_rnf      ! river runoff 
    23    USE sbcrnf , ONLY : sbc_rnf_div ! river runoff  
     23   USE sbcrnf , ONLY : sbc_rnf_div ! river runoff 
    2424   USE isf_oce, ONLY : ln_isf      ! ice shelf 
    2525   USE isfhdiv, ONLY : isf_hdiv    ! ice shelf 
    26 #if defined key_asminc    
     26#if defined key_asminc 
    2727   USE asminc          ! Assimilation increment 
    2828#endif 
    … …  
    4040   !! * Substitutions 
    4141#  include "do_loop_substitute.h90" 
     42#  include "domzgr_substitute.h90" 
    4243   !!---------------------------------------------------------------------- 
    4344   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    44    !! $Id$  
     45   !! $Id$ 
    4546   !! Software governed by the CeCILL license (see ./LICENSE) 
    4647   !!---------------------------------------------------------------------- 
    … …  
    5051      !!---------------------------------------------------------------------- 
    5152      !!                  ***  ROUTINE div_hor  *** 
    52       !!                     
     53      !! 
    5354      !! ** Purpose :   compute the horizontal divergence at now time-step 
    5455      !! 
    5556      !! ** Method  :   the now divergence is computed as : 
    5657      !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    57       !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla)  
     58      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 
    5859      !! 
    5960      !! ** Action  : - update hdiv, the now horizontal divergence 
    … …  
    7879      DO_3D_00_00( 1, jpkm1 ) 
    7980         hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      & 
    80             &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      & 
    81             &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      & 
    82             &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   & 
     81            &              - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      & 
     82            &              + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      & 
     83            &              - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   & 
    8384            &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    8485      END_3D 
    … …  
    9596      IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )                     !==  runoffs    ==!   (update hdiv field) 
    9697      ! 
    97 #if defined key_asminc  
     98#if defined key_asminc 
    9899      IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, Kbb, Kmm, hdiv )   !==  SSH assimilation  ==!   (update hdiv field) 
    99       !  
     100      ! 
    100101#endif 
    101102      ! 
    … …  
    107108      ! 
    108109   END SUBROUTINE div_hor 
    109     
     110 
    110111   !!====================================================================== 
    111112END MODULE divhor 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynvor.F90

    r12377 r12590  
    1515   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme 
    1616   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase 
    17    !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity  
     17   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity 
    1818   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory 
    1919   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 
    … …  
    7070   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme 
    7171 
    72    INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
     72   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
    7373   !                               ! associated indices: 
    7474   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
    … …  
    7979 
    8080   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation 
    81    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -  
     81   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       - 
    8282   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation 
    83    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -  
    84     
     83   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       - 
     84 
    8585   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
    8686   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8 
    8787   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12 
    88     
     88 
    8989   !! * Substitutions 
    9090#  include "do_loop_substitute.h90" 
     91#  include "domzgr_substitute.h90" 
     92 
    9193   !!---------------------------------------------------------------------- 
    9294   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    103105      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 
    104106      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    105       !!               and planetary vorticity trends) and send them to trd_dyn  
     107      !!               and planetary vorticity trends) and send them to trd_dyn 
    106108      !!               for futher diagnostics (l_trddyn=T) 
    107109      !!---------------------------------------------------------------------- 
    … …  
    193195      !!                  ***  ROUTINE vor_enT  *** 
    194196      !! 
    195       !! ** Purpose :   Compute the now total vorticity trend and add it to  
     197      !! ** Purpose :   Compute the now total vorticity trend and add it to 
    196198      !!      the general trend of the momentum equation. 
    197199      !! 
    198       !! ** Method  :   Trend evaluated using now fields (centered in time)  
     200      !! ** Method  :   Trend evaluated using now fields (centered in time) 
    199201      !!       and t-point evaluation of vorticity (planetary and relative). 
    200202      !!       conserves the horizontal kinetic energy. 
    201       !!         The general trend of momentum is increased due to the vorticity  
     203      !!         The general trend of momentum is increased due to the vorticity 
    202204      !!       term which is given by: 
    203205      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ] 
    … …  
    233235                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    234236            END_2D 
    235             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
     237            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
    236238               DO_2D_10_10 
    237239                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    … …  
    248250                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
    249251            END_2D 
    250             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
     252            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
    251253               DO_2D_10_10 
    252254                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    … …  
    269271            DO_2D_01_01 
    270272               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
    271                   &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
     273                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) & 
     274                  &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    272275            END_2D 
    273276         CASE ( np_MET )                           !* metric term 
    274277            DO_2D_01_01 
    275                zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
    276                   &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t(ji,jj,jk,Kmm) 
     278               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)     & 
     279                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) & 
     280                  &             * e3t(ji,jj,jk,Kmm) 
    277281            END_2D 
    278282         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    279283            DO_2D_01_01 
    280                zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
    281                   &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
     284               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)      & 
     285                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) & 
     286                  &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    282287            END_2D 
    283288         CASE ( np_CME )                           !* Coriolis + metric 
    284289            DO_2D_01_01 
    285                zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
    286                     &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
    287                     &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t(ji,jj,jk,Kmm) 
     290               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                             & 
     291                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)    & 
     292                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) & 
     293                    &          * e3t(ji,jj,jk,Kmm) 
    288294            END_2D 
    289295         CASE DEFAULT                                             ! error 
    … …  
    298304               ! 
    299305            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    & 
    300                &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &  
    301                &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )  
     306               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
     307               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
    302308         END_2D 
    303309         !                                             ! =============== 
    … …  
    311317      !!                  ***  ROUTINE vor_ene  *** 
    312318      !! 
    313       !! ** Purpose :   Compute the now total vorticity trend and add it to  
     319      !! ** Purpose :   Compute the now total vorticity trend and add it to 
    314320      !!      the general trend of the momentum equation. 
    315321      !! 
    316       !! ** Method  :   Trend evaluated using now fields (centered in time)  
     322      !! ** Method  :   Trend evaluated using now fields (centered in time) 
    317323      !!       and the Sadourny (1975) flux form formulation : conserves the 
    318324      !!       horizontal kinetic energy. 
    319       !!         The general trend of momentum is increased due to the vorticity  
     325      !!         The general trend of momentum is increased due to the vorticity 
    320326      !!       term which is given by: 
    321327      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ] 
    … …  
    350356         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
    351357         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    352             zwz(:,:) = ff_f(:,:)  
     358            zwz(:,:) = ff_f(:,:) 
    353359         CASE ( np_RVO )                           !* relative vorticity 
    354360            DO_2D_10_10 
    … …  
    396402            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    397403            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    398             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     404            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    399405         END_2D 
    400406         !                                             ! =============== 
    … …  
    446452         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
    447453         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    448             zwz(:,:) = ff_f(:,:)  
     454            zwz(:,:) = ff_f(:,:) 
    449455         CASE ( np_RVO )                           !* relative vorticity 
    450456            DO_2D_10_10 
    … …  
    504510      !!                ***  ROUTINE vor_een  *** 
    505511      !! 
    506       !! ** Purpose :   Compute the now total vorticity trend and add it to  
     512      !! ** Purpose :   Compute the now total vorticity trend and add it to 
    507513      !!      the general trend of the momentum equation. 
    508514      !! 
    509       !! ** Method  :   Trend evaluated using now fields (centered in time)  
    510       !!      and the Arakawa and Lamb (1980) flux form formulation : conserves  
     515      !! ** Method  :   Trend evaluated using now fields (centered in time) 
     516      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves 
    511517      !!      both the horizontal kinetic energy and the potential enstrophy 
    512518      !!      when horizontal divergence is zero (see the NEMO documentation) 
    … …  
    545551         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    546552            DO_2D_10_10 
    547                ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
    548                   &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     553               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     554                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     555                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     556                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
    549557               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f 
    550558               ELSE                       ;   z1_e3f(ji,jj) = 0._wp 
    … …  
    553561         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
    554562            DO_2D_10_10 
    555                ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
    556                   &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     563               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     564                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     565                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     566                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
    557567               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
    558568                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  ) 
    … …  
    644654      !!                ***  ROUTINE vor_eeT  *** 
    645655      !! 
    646       !! ** Purpose :   Compute the now total vorticity trend and add it to  
     656      !! ** Purpose :   Compute the now total vorticity trend and add it to 
    647657      !!      the general trend of the momentum equation. 
    648658      !! 
    649       !! ** Method  :   Trend evaluated using now fields (centered in time)  
    650       !!      and the Arakawa and Lamb (1980) vector form formulation using  
     659      !! ** Method  :   Trend evaluated using now fields (centered in time) 
     660      !!      and the Arakawa and Lamb (1980) vector form formulation using 
    651661      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een). 
    652       !!      The change consists in  
     662      !!      The change consists in 
    653663      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs). 
    654664      !! 
    … …  
    667677      REAL(wp) ::   zua, zva       ! local scalars 
    668678      REAL(wp) ::   zmsk, z1_e3t   ! local scalars 
    669       REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy  
     679      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy 
    670680      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse 
    671681      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz 
    … …  
    827837      ! 
    828838      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 
    829       !                       
     839      ! 
    830840      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot) 
    831841      ncor = np_COR                       ! planetary vorticity 
    … …  
    836846         ntot = np_COR        !     -         - 
    837847      CASE( np_VEC_c2  ) 
    838          IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity'  
     848         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
    839849         nrvm = np_RVO        ! relative vorticity 
    840          ntot = np_CRV        ! relative + planetary vorticity          
     850         ntot = np_CRV        ! relative + planetary vorticity 
    841851      CASE( np_FLX_c2 , np_FLX_ubs  ) 
    842852         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term' 
    … …  
    863873         ! 
    864874      END SELECT 
    865        
     875 
    866876      IF(lwp) THEN                   ! Print the choice 
    867877         WRITE(numout,*) 
    … …  
    873883         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' 
    874884         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)' 
    875          END SELECT          
     885         END SELECT 
    876886      ENDIF 
    877887      ! 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90

    r12377 r12590  
    1 MODULE sshwzv    
     1MODULE sshwzv 
    22   !!============================================================================== 
    33   !!                       ***  MODULE  sshwzv  *** 
    … …  
    55   !!============================================================================== 
    66   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code 
    7    !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA  
     7   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA 
    88   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    99   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    … …  
    2020   USE oce            ! ocean dynamics and tracers variables 
    2121   USE isf_oce        ! ice shelf 
    22    USE dom_oce        ! ocean space and time domain variables  
     22   USE dom_oce        ! ocean space and time domain variables 
    2323   USE sbc_oce        ! surface boundary condition: ocean 
    2424   USE domvvl         ! Variable volume 
    … …  
    3131#endif 
    3232   ! 
    33    USE iom  
     33   USE iom 
    3434   USE in_out_manager ! I/O manager 
    3535   USE restart        ! only for lrst_oce 
    … …  
    5050   !! * Substitutions 
    5151#  include "do_loop_substitute.h90" 
     52#  include "domzgr_substitute.h90" 
     53 
    5254   !!---------------------------------------------------------------------- 
    5355   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    6062      !!---------------------------------------------------------------------- 
    6163      !!                ***  ROUTINE ssh_nxt  *** 
    62       !!                    
     64      !! 
    6365      !! ** Purpose :   compute the after ssh (ssh(Kaa)) 
    6466      !! 
    … …  
    7476      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index 
    7577      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height 
    76       !  
     78      ! 
    7779      INTEGER  ::   jk            ! dummy loop indice 
    7880      REAL(wp) ::   z2dt, zcoef   ! local scalars 
    … …  
    108110      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to 
    109111      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    110       !  
     112      ! 
    111113      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    112114      ! 
    … …  
    131133   END SUBROUTINE ssh_nxt 
    132134 
    133     
     135 
    134136   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa ) 
    135137      !!---------------------------------------------------------------------- 
    136138      !!                ***  ROUTINE wzv  *** 
    137       !!                    
     139      !! 
    138140      !! ** Purpose :   compute the now vertical velocity 
    139141      !! 
    140       !! ** Method  : - Using the incompressibility hypothesis, the vertical  
    141       !!      velocity is computed by integrating the horizontal divergence   
     142      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
     143      !!      velocity is computed by integrating the horizontal divergence 
    142144      !!      from the bottom to the surface minus the scale factor evolution. 
    143145      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
    … …  
    172174      ! 
    173175      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
    174          ALLOCATE( zhdiv(jpi,jpj,jpk) )  
     176         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
    175177         ! 
    176178         DO jk = 1, jpkm1 
    … …  
    186188         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    187189            ! computation of w 
    188             pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk)    & 
    189                &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )     ) * tmask(:,:,jk) 
     190            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    & 
     191               &                         + zhdiv(:,:,jk)                       & 
     192               &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          & 
     193               &                                    - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk) 
    190194         END DO 
    191195         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 
    192          DEALLOCATE( zhdiv )  
     196         DEALLOCATE( zhdiv ) 
    193197      ELSE   ! z_star and linear free surface cases 
    194198         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    195199            ! computation of w 
    196             pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 & 
    197                &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk) 
     200            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    & 
     201               &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          & 
     202               &                                    - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk) 
    198203         END DO 
    199204      ENDIF 
    … …  
    205210      ENDIF 
    206211      ! 
    207 #if defined key_agrif  
    208       IF( .NOT. AGRIF_Root() ) THEN  
    209          IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east  
    210          IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west  
    211          IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north  
    212          IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south  
    213       ENDIF  
    214 #endif  
     212#if defined key_agrif 
     213      IF( .NOT. AGRIF_Root() ) THEN 
     214         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east 
     215         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west 
     216         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north 
     217         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south 
     218      ENDIF 
     219#endif 
    215220      ! 
    216221      IF( ln_timing )   CALL timing_stop('wzv') 
    … …  
    274279      !!---------------------------------------------------------------------- 
    275280      !!                ***  ROUTINE wAimp  *** 
    276       !!                    
     281      !! 
    277282      !! ** Purpose :   compute the Courant number and partition vertical velocity 
    278283      !!                if a proportion needs to be treated implicitly 
    279284      !! 
    280       !! ** Method  : -  
     285      !! ** Method  : - 
    281286      !! 
    282287      !! ** action  :   ww      : now vertical velocity (to be handled explicitly) 
    … …  
    284289      !! 
    285290      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent 
    286       !!              implicit scheme for vertical advection in oceanic modeling.  
     291      !!              implicit scheme for vertical advection in oceanic modeling. 
    287292      !!              Ocean Modelling, 91, 38-69. 
    288293      !!---------------------------------------------------------------------- 
    … …  
    312317            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
    313318            ! 2*rdt and not r2dt (for restartability) 
    314             Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )                       &   
     319            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )                       & 
    315320               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   & 
    316321               &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) )   & 
    … …  
    325330            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
    326331            ! 2*rdt and not r2dt (for restartability) 
    327             Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )   &  
     332            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )   & 
    328333               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   & 
    329334               &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   & 
    … …  
    344349            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) 
    345350! alt: 
    346 !                  IF ( ww(ji,jj,jk) > 0._wp ) THEN  
    347 !                     zCu =  Cu_adv(ji,jj,jk)  
     351!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN 
     352!                     zCu =  Cu_adv(ji,jj,jk) 
    348353!                  ELSE 
    349354!                     zCu =  Cu_adv(ji,jj,jk-1) 
    350 !                  ENDIF  
     355!                  ENDIF 
    351356            ! 
    352357            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit 
    … …  
    365370            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl 
    366371         END_3D 
    367          Cu_adv(:,:,1) = 0._wp  
     372         Cu_adv(:,:,1) = 0._wp 
    368373      ELSE 
    369374         ! Fully explicit everywhere 
    … …  
    371376         wi    (:,:,:) = 0._wp 
    372377      ENDIF 
    373       CALL iom_put("wimp",wi)  
     378      CALL iom_put("wimp",wi) 
    374379      CALL iom_put("wi_cff",Cu_adv) 
    375380      CALL iom_put("wexp",ww) 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/FLO/floblk.F90

    r12377 r12590  
    2020   PUBLIC   flo_blk    ! routine called by floats.F90 
    2121 
     22#  include "domzgr_substitute.h90" 
     23 
    2224   !!---------------------------------------------------------------------- 
    2325   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    24    !! $Id$  
     26   !! $Id$ 
    2527   !! Software governed by the CeCILL license (see ./LICENSE) 
    2628   !!---------------------------------------------------------------------- 
    … …  
    3032      !!--------------------------------------------------------------------- 
    3133      !!                  ***  ROUTINE flo_blk  *** 
    32       !!            
     34      !! 
    3335      !! ** Purpose :   Compute the geographical position,latitude, longitude 
    3436      !!      and depth of each float at each time step. 
    35       !!  
     37      !! 
    3638      !! ** Method  :   The position of a float is computed with Bruno Blanke 
    3739      !!      algorithm. We need to know the velocity field, the old positions 
    … …  
    4749         zuoutfl,zvoutfl,zwoutfl,   &     ! transport across the ouput face 
    4850         zvol,                      &     ! volume of the mesh 
    49          zsurfz,                    &     ! surface of the face of the mesh  
     51         zsurfz,                    &     ! surface of the face of the mesh 
    5052         zind 
    5153 
    … …  
    5355 
    5456      INTEGER  , DIMENSION ( jpnfl )  ::   iil, ijl, ikl                   ! index of nearest mesh 
    55       INTEGER  , DIMENSION ( jpnfl )  ::   iiloc , ijloc               
     57      INTEGER  , DIMENSION ( jpnfl )  ::   iiloc , ijloc 
    5658      INTEGER  , DIMENSION ( jpnfl )  ::   iiinfl, ijinfl, ikinfl          ! index of input mesh of the float. 
    5759      INTEGER  , DIMENSION ( jpnfl )  ::   iioutfl, ijoutfl, ikoutfl       ! index of output mesh of the float. 
    58       REAL(wp) , DIMENSION ( jpnfl )  ::   zgifl, zgjfl, zgkfl             ! position of floats, index on  
     60      REAL(wp) , DIMENSION ( jpnfl )  ::   zgifl, zgjfl, zgkfl             ! position of floats, index on 
    5961      !                                                                         ! velocity mesh. 
    6062      REAL(wp) , DIMENSION ( jpnfl )  ::    ztxfl, ztyfl, ztzfl            ! time for a float to quit the mesh 
    61       !                                                                         ! across one of the face x,y and z  
    62       REAL(wp) , DIMENSION ( jpnfl )  ::    zttfl                          ! time for a float to quit the mesh  
    63       REAL(wp) , DIMENSION ( jpnfl )  ::    zagefl                         ! time during which, trajectorie of  
     63      !                                                                         ! across one of the face x,y and z 
     64      REAL(wp) , DIMENSION ( jpnfl )  ::    zttfl                          ! time for a float to quit the mesh 
     65      REAL(wp) , DIMENSION ( jpnfl )  ::    zagefl                         ! time during which, trajectorie of 
    6466      !                                                                         ! the float has been computed 
    65       REAL(wp) , DIMENSION ( jpnfl )  ::   zagenewfl                       ! new age of float after calculation  
     67      REAL(wp) , DIMENSION ( jpnfl )  ::   zagenewfl                       ! new age of float after calculation 
    6668      !                                                                         ! of new position 
    6769      REAL(wp) , DIMENSION ( jpnfl )  ::   zufl, zvfl, zwfl                ! interpolated vel. at float position 
    … …  
    7779 
    7880      ! Initialisation of parameters 
    79        
     81 
    8082      DO jfl = 1, jpnfl 
    8183         ! ages of floats are put at zero 
    8284         zagefl(jfl) = 0. 
    83          ! index on the velocity grid  
    84          ! We considere k coordinate negative, with this transformation  
    85          ! the computation in the 3 direction is the same.  
     85         ! index on the velocity grid 
     86         ! We considere k coordinate negative, with this transformation 
     87         ! the computation in the 3 direction is the same. 
    8688         zgifl(jfl) = tpifl(jfl) - 0.5 
    8789         zgjfl(jfl) = tpjfl(jfl) - 0.5 
    8890         zgkfl(jfl) = MIN(-1.,-(tpkfl(jfl))) 
    89          ! surface drift every 10 days  
     91         ! surface drift every 10 days 
    9092         IF( ln_argo ) THEN 
    9193            IF( MOD(kt,150) >= 146 .OR. MOD(kt,150) == 0 )  zgkfl(jfl) = -1. 
    … …  
    9698         ikl(jfl) =     INT(zgkfl(jfl)) 
    9799      END DO 
    98         
     100 
    99101      iloop = 0 
    100102222   DO jfl = 1, jpnfl 
    … …  
    104106            iiloc(jfl) = iil(jfl) - mig(1) + 1 
    105107            ijloc(jfl) = ijl(jfl) - mjg(1) + 1 
    106 # else  
     108# else 
    107109            iiloc(jfl) = iil(jfl) 
    108110            ijloc(jfl) = ijl(jfl) 
    109111# endif 
    110              
    111             ! compute the transport across the mesh where the float is.             
    112 !!bug (gm) change e3t into e3. but never checked  
    113             zsurfx(1) = e2u(iiloc(jfl)-1,ijloc(jfl)  ) * e3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl),Kmm) 
    114             zsurfx(2) = e2u(iiloc(jfl)  ,ijloc(jfl)  ) * e3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm) 
    115             zsurfy(1) = e1v(iiloc(jfl)  ,ijloc(jfl)-1) * e3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl),Kmm) 
    116             zsurfy(2) = e1v(iiloc(jfl)  ,ijloc(jfl)  ) * e3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm) 
     112 
     113            ! compute the transport across the mesh where the float is. 
     114!!bug (gm) change e3t into e3. but never checked 
     115            zsurfx(1) = e2u(iiloc(jfl)-1,ijloc(jfl)  )    & 
     116            &         * e3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl),Kmm) 
     117            zsurfx(2) = e2u(iiloc(jfl)  ,ijloc(jfl)  )    & 
     118            &         * e3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm) 
     119            zsurfy(1) = e1v(iiloc(jfl)  ,ijloc(jfl)-1)    & 
     120            &         * e3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl),Kmm) 
     121            zsurfy(2) = e1v(iiloc(jfl)  ,ijloc(jfl)  )    & 
     122            &         * e3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm) 
    117123 
    118124            ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too. 
    … …  
    129135            zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl)   )   & 
    130136               &   +  ww(iiloc(jfl),ijloc(jfl),- ikl(jfl)   ) )/2. *  zsurfz*nisobfl(jfl) 
    131              
    132             ! interpolation of velocity field on the float initial position             
     137 
     138            ! interpolation of velocity field on the float initial position 
    133139            zufl(jfl)=  zuinfl  + ( zgifl(jfl) - float(iil(jfl)-1) ) * ( zuoutfl - zuinfl) 
    134140            zvfl(jfl)=  zvinfl  + ( zgjfl(jfl) - float(ijl(jfl)-1) ) * ( zvoutfl - zvinfl) 
    135141            zwfl(jfl)=  zwinfl  + ( zgkfl(jfl) - float(ikl(jfl)-1) ) * ( zwoutfl - zwinfl) 
    136              
     142 
    137143            ! faces of input and output 
    138144            ! u-direction 
    … …  
    147153               iiinfl (jfl) = iil(jfl) - 1 
    148154            ENDIF 
    149             ! v-direction        
     155            ! v-direction 
    150156            IF( zvfl(jfl) < 0. ) THEN 
    151157               ijoutfl(jfl) = ijl(jfl) - 1. 
    … …  
    169175               ikinfl (jfl) = ikl(jfl) - 1. 
    170176            ENDIF 
    171              
     177 
    172178            ! compute the time to go out the mesh across a face 
    173179            ! u-direction 
    … …  
    203209               ENDIF 
    204210            ENDIF 
    205             ! w-direction         
    206             IF( nisobfl(jfl) == 1. ) THEN  
     211            ! w-direction 
     212            IF( nisobfl(jfl) == 1. ) THEN 
    207213               zwdfl (jfl) = zwoutfl - zwinfl 
    208214               zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl)) 
    … …  
    221227               ENDIF 
    222228            ENDIF 
    223              
     229 
    224230            ! the time to go leave the mesh is the smallest time 
    225                     
    226             IF( nisobfl(jfl) == 1. ) THEN  
     231 
     232            IF( nisobfl(jfl) == 1. ) THEN 
    227233               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl)) 
    228234            ELSE 
    … …  
    231237            ! new age of the FLOAT 
    232238            zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol 
    233             ! test to know if the "age" of the float is not bigger than the  
     239            ! test to know if the "age" of the float is not bigger than the 
    234240            ! time step 
    235241            IF( zagenewfl(jfl) > rdt ) THEN 
    … …  
    237243               zagenewfl(jfl) = rdt 
    238244            ENDIF 
    239              
     245 
    240246            ! In the "minimal" direction we compute the index of new mesh 
    241247            ! on i-direction 
    … …  
    250256               iiinfl(jfl) = ind 
    251257            ELSE 
    252                IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN  
     258               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN 
    253259                  zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl)    & 
    254260                     &       * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) /  zudfl(jfl) 
    … …  
    268274               ijinfl(jfl) = ind 
    269275            ELSE 
    270                IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN  
     276               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN 
    271277                  zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl)   & 
    272278                     &       * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) /  zvdfl(jfl) 
    … …  
    287293                  ikinfl(jfl) = ind 
    288294               ELSE 
    289                   IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN  
     295                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN 
    290296                     zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl)    & 
    291297                        &       * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) /  zwdfl(jfl) 
    … …  
    295301               ENDIF 
    296302            ENDIF 
    297              
     303 
    298304            ! coordinate of the new point on the temperature grid 
    299              
     305 
    300306            iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl)) 
    301307            ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl)) 
    … …  
    306312!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl) 
    307313!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl), 
    308 !!Alexcadm     .  zgjfl(jfl)  
     314!!Alexcadm     .  zgjfl(jfl) 
    309315!!Alexcadm  IF (jfl == 910) write(*,*)'Flotteur 910', 
    310316!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl) 
    … …  
    312318!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl) 
    313319!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl), 
    314 !!Alexcadm     .  zgjfl(jfl)  
     320!!Alexcadm     .  zgjfl(jfl) 
    315321            ! reinitialisation of the age of FLOAT 
    316322            zagefl(jfl) = zagenewfl(jfl) 
    … …  
    327333# endif 
    328334      END DO 
    329        
     335 
    330336      ! synchronisation 
    331337      CALL mpp_sum( 'floblk', zgifl , jpnfl )   ! sums over the global domain 
    … …  
    335341      CALL mpp_sum( 'floblk', iil   , jpnfl ) 
    336342      CALL mpp_sum( 'floblk', ijl   , jpnfl ) 
    337        
     343 
    338344      ! Test to know if a  float hasn't integrated enought time 
    339345      IF( ln_argo ) THEN 
    … …  
    361367!!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl) 
    362368      IF( ifin == 0 ) THEN 
    363          iloop = iloop + 1  
     369         iloop = iloop + 1 
    364370         GO TO 222 
    365371      ENDIF 
    … …  
    369375 
    370376   !!====================================================================== 
    371 END MODULE floblk  
     377END MODULE floblk 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcrnf.F90

    r12377 r12590  
    3434   PUBLIC   sbc_rnf_alloc ! called in sbcmod module 
    3535   PUBLIC   sbc_rnf_init  ! called in sbcmod module 
    36     
     36 
    3737   !                                                !!* namsbc_rnf namelist * 
    3838   CHARACTER(len=100)         ::   cn_dir            !: Root directory for location of rnf files 
    … …  
    5858   LOGICAL , PUBLIC ::   l_rnfcpl = .false.   !: runoffs recieved from oasis 
    5959   INTEGER , PUBLIC ::   nkrnf = 0            !: nb of levels over which Kz is increased at river mouths 
    60     
     60 
    6161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rnfmsk              !: river mouth mask (hori.) 
    6262   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   rnfmsk_z            !: river mouth mask (vert.) 
    6363   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   h_rnf               !: depth of runoff in m 
    6464   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nk_rnf              !: depth of runoff in model levels 
    65    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rnf_tsc_b, rnf_tsc  !: before and now T & S runoff contents   [K.m/s & PSU.m/s]    
     65   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rnf_tsc_b, rnf_tsc  !: before and now T & S runoff contents   [K.m/s & PSU.m/s] 
    6666 
    6767   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_rnf       ! structure: river runoff (file information, fields read) 
    6868   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_i_rnf     ! structure: iceberg flux (file information, fields read) 
    69    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read)   
    70    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)   
    71   
     69   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read) 
     70   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read) 
     71 
    7272   !! * Substitutions 
    7373#  include "do_loop_substitute.h90" 
     74#  include "domzgr_substitute.h90" 
    7475   !!---------------------------------------------------------------------- 
    7576   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    227228      ELSE                       !==   runoff put only at the surface   ==! 
    228229         h_rnf (:,:)   = e3t (:,:,1,Kmm)        ! update h_rnf to be depth of top box 
    229          phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / e3t(:,:,1,Kmm) 
     230         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:)+rnf_b(:,:) ) * zfact * r1_rau0 / e3t(:,:,1,Kmm) 
    230231      ENDIF 
    231232      ! 
    … …  
    249250      INTEGER           ::   ios           ! Local integer output status for namelist read 
    250251      INTEGER           ::   nbrec         ! temporary integer 
    251       REAL(wp)          ::   zacoef   
    252       REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl     
     252      REAL(wp)          ::   zacoef 
     253      REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl 
    253254      !! 
    254255      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb,   & 
    … …  
    261262      IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_rnf_alloc : unable to allocate arrays' ) 
    262263      ! 
    263       IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths  
     264      IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths 
    264265         ln_rnf_mouth  = .FALSE.                   ! default definition needed for example by sbc_ssr or by tra_adv_muscl 
    265266         nkrnf         = 0 
    … …  
    297298      !                                   ! ================== 
    298299      ! 
    299       IF( .NOT. l_rnfcpl ) THEN                     
     300      IF( .NOT. l_rnfcpl ) THEN 
    300301         ALLOCATE( sf_rnf(1), STAT=ierror )         ! Create sf_rnf structure (runoff inflow) 
    301302         IF(lwp) WRITE(numout,*) 
    … …  
    352353         IF(lwp) WRITE(numout,*) '   ==>>>   runoffs depth read in a file' 
    353354         rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname ) 
    354          IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year  
    355             IF( sn_dep_rnf%cltype == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month  
     355         IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year 
     356            IF( sn_dep_rnf%cltype == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month 
    356357         ENDIF 
    357358         CALL iom_open ( rn_dep_file, inum )                           ! open file 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcssm.F90

    r12377 r12590  
    1010 
    1111   !!---------------------------------------------------------------------- 
    12    !!   sbc_ssm       : calculate sea surface mean currents, temperature,   
     12   !!   sbc_ssm       : calculate sea surface mean currents, temperature, 
    1313   !!                   and salinity over nn_fsbc time-step 
    1414   !!---------------------------------------------------------------------- 
    … …  
    3131 
    3232   LOGICAL, SAVE ::   l_ssm_mean = .FALSE.   ! keep track of whether means have been read from restart file 
    33     
     33 
     34#  include "domzgr_substitute.h90" 
    3435   !!---------------------------------------------------------------------- 
    3536   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    4243      !!--------------------------------------------------------------------- 
    4344      !!                     ***  ROUTINE sbc_oce  *** 
    44       !!                      
     45      !! 
    4546      !! ** Purpose :   provide ocean surface variable to sea-surface boundary 
    46       !!                condition computation  
    47       !!                 
    48       !! ** Method  :   compute mean surface velocity (2 components at U and  
     47      !!                condition computation 
     48      !! 
     49      !! ** Method  :   compute mean surface velocity (2 components at U and 
    4950      !!      V-points) [m/s], temperature [Celsius] and salinity [psu] over 
    5051      !!      the periode (kt - nn_fsbc) to kt 
    … …  
    200201         ! 
    201202      ELSE 
    202          !                
     203         ! 
    203204         IF(lwp) WRITE(numout,*) 
    204205         IF(lwp) WRITE(numout,*) 'sbc_ssm_init : sea surface mean fields' 
    … …  
    222223            ! 
    223224            IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs 
    224                IF(lwp) WRITE(numout,*) '   restart with a change in the frequency of mean from ', zf_sbc, ' to ', nn_fsbc  
    225                zcoef = REAL( nn_fsbc - 1, wp ) / zf_sbc  
    226                ssu_m(:,:) = zcoef * ssu_m(:,:)  
     225               IF(lwp) WRITE(numout,*) '   restart with a change in the frequency of mean from ', zf_sbc, ' to ', nn_fsbc 
     226               zcoef = REAL( nn_fsbc - 1, wp ) / zf_sbc 
     227               ssu_m(:,:) = zcoef * ssu_m(:,:) 
    227228               ssv_m(:,:) = zcoef * ssv_m(:,:) 
    228229               sst_m(:,:) = zcoef * sst_m(:,:) 
    … …  
    252253      ENDIF 
    253254      ! 
    254       IF( .NOT. ln_traqsr )   fraqsr_1lev(:,:) = 1._wp   ! default definition: qsr 100% in the fisrt level  
     255      IF( .NOT. ln_traqsr )   fraqsr_1lev(:,:) = 1._wp   ! default definition: qsr 100% in the fisrt level 
    255256      ! 
    256257      IF( lwxios.AND.nn_fsbc > 1 ) THEN 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_cen.F90

    r12377 r12590  
    1313   USE dom_oce        ! ocean space and time domain 
    1414   USE eosbn2         ! equation of state 
    15    USE traadv_fct     ! acces to routine interp_4th_cpt  
     15   USE traadv_fct     ! acces to routine interp_4th_cpt 
    1616   USE trd_oce        ! trends: ocean variables 
    17    USE trdtra         ! trends manager: tracers  
     17   USE trdtra         ! trends manager: tracers 
    1818   USE diaptr         ! poleward transport diagnostics 
    1919   USE diaar5         ! AR5 diagnostics 
    … …  
    2828 
    2929   PUBLIC   tra_adv_cen   ! called by traadv.F90 
    30     
     30 
    3131   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6 
    3232 
    … …  
    3737   !! * Substitutions 
    3838#  include "do_loop_substitute.h90" 
     39#  include "domzgr_substitute.h90" 
    3940   !!---------------------------------------------------------------------- 
    4041   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    4546 
    4647   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW,     & 
    47       &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )  
     48      &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 
    4849      !!---------------------------------------------------------------------- 
    4950      !!                  ***  ROUTINE tra_adv_cen  *** 
    50       !!                  
     51      !! 
    5152      !! ** Purpose :   Compute the now trend due to the advection of tracers 
    5253      !!      and add it to the general trend of passive tracer equations. 
    5354      !! 
    5455      !! ** Method  :   The advection is evaluated by a 2nd or 4th order scheme 
    55       !!               using now fields (leap-frog scheme).  
     56      !!               using now fields (leap-frog scheme). 
    5657      !!       kn_cen_h = 2  ==>> 2nd order centered scheme on the horizontal 
    5758      !!                = 4  ==>> 4th order    -        -       -      - 
    … …  
    9091      l_ptr = .FALSE. 
    9192      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )       l_trd = .TRUE. 
    92       IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )    l_ptr = .TRUE.  
     93      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )    l_ptr = .TRUE. 
    9394      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    9495         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
    9596      ! 
    96       !                     
     97      ! 
    9798      zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers 
    9899      zwz(:,:,jpk) = 0._wp 
    … …  
    150151            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean) 
    151152               DO_2D_11_11 
    152                   zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)  
     153                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
    153154               END_2D 
    154155            ELSE                                   ! no ice-shelf cavities (only ocean surface) 
    … …  
    156157            ENDIF 
    157158         ENDIF 
    158          !                
     159         ! 
    159160         DO_3D_00_00( 1, jpkm1 ) 
    160161            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)    & 
    161162               &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    & 
    162163               &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    & 
    163                &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     164               &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) & 
     165               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    164166         END_3D 
    165167         !                             ! trend diagnostics 
    … …  
    169171            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) ) 
    170172         END IF 
    171          !                                 ! "Poleward" heat and salt transports  
     173         !                                 ! "Poleward" heat and salt transports 
    172174         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 
    173175         !                                 !  heat and salt transport 
    … …  
    177179      ! 
    178180   END SUBROUTINE tra_adv_cen 
    179     
     181 
    180182   !!====================================================================== 
    181183END MODULE traadv_cen 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_fct.F90

    r12377 r12590  
    1010   !!  tra_adv_fct    : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme 
    1111   !!                   with sub-time-stepping in the vertical direction 
    12    !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     12   !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm 
    1313   !!  interp_4th_cpt : 4th order compact scheme for the vertical component of the advection 
    1414   !!---------------------------------------------------------------------- 
    … …  
    2424   ! 
    2525   USE in_out_manager ! I/O manager 
    26    USE iom            !  
     26   USE iom            ! 
    2727   USE lib_mpp        ! MPP library 
    28    USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    29    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     28   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
     29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3030 
    3131   IMPLICIT NONE 
    … …  
    4646   !! * Substitutions 
    4747#  include "do_loop_substitute.h90" 
     48#  include "domzgr_substitute.h90" 
    4849   !!---------------------------------------------------------------------- 
    4950   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5758      !!---------------------------------------------------------------------- 
    5859      !!                  ***  ROUTINE tra_adv_fct  *** 
    59       !!  
     60      !! 
    6061      !! **  Purpose :   Compute the now trend due to total advection of tracers 
    6162      !!               and add it to the general trend of tracer equations 
    … …  
    6364      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction 
    6465      !!               (choice through the value of kn_fct) 
    65       !!               - on the vertical the 4th order is a compact scheme  
    66       !!               - corrected flux (monotonic correction)  
     66      !!               - on the vertical the 4th order is a compact scheme 
     67      !!               - corrected flux (monotonic correction) 
    6768      !! 
    6869      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
    … …  
    8182      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    8283      ! 
    83       INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices   
     84      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices 
    8485      REAL(wp) ::   ztra                                     ! local scalar 
    8586      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      - 
    … …  
    102103      ll_zAimp = .FALSE. 
    103104      IF( ( cdtype == 'TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
    104       IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE.  
     105      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE. 
    105106      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  & 
    106107         &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
    … …  
    111112      ENDIF 
    112113      ! 
    113       IF( l_ptr ) THEN   
     114      IF( l_ptr ) THEN 
    114115         ALLOCATE( zptry(jpi,jpj,jpk) ) 
    115116         zptry(:,:,:) = 0._wp 
    116117      ENDIF 
    117118      !                          ! surface & bottom value : flux set to zero one for all 
    118       zwz(:,:, 1 ) = 0._wp             
     119      zwz(:,:, 1 ) = 0._wp 
    119120      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp 
    120121      ! 
    121       zwi(:,:,:) = 0._wp         
     122      zwi(:,:,:) = 0._wp 
    122123      ! 
    123124      ! If adaptive vertical advection, check if it is needed on this PE at this time 
    … …  
    129130         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 
    130131         DO_3D_00_00( 1, jpkm1 ) 
    131             zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs) 
     132            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )   & 
     133            &                               / e3t(ji,jj,jk,Krhs) 
    132134            zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs) 
    133135            zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) 
    … …  
    138140         ! 
    139141         !        !==  upstream advection with initial mass fluxes & intermediate update  ==! 
    140          !                    !* upstream tracer flux in the i and j direction  
     142         !                    !* upstream tracer flux in the i and j direction 
    141143         DO_3D_10_10( 1, jpkm1 ) 
    142144            ! upstream scheme 
    … …  
    157159            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface 
    158160               DO_2D_11_11 
    159                   zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface  
     161                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
    160162               END_2D 
    161163            ELSE                             ! no cavities: only at the ocean surface 
    … …  
    163165            ENDIF 
    164166         ENDIF 
    165          !                
     167         ! 
    166168         DO_3D_00_00( 1, jpkm1 ) 
    167169            !                             ! total intermediate advective trends 
    … …  
    170172               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
    171173            !                             ! update and guess with monotonic sheme 
    172             pt(ji,jj,jk,jn,Krhs) =                     pt(ji,jj,jk,jn,Krhs) +        ztra   / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
    173             zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     174            pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   & 
     175               &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 
     176            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 
     177               &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
    174178         END_3D 
    175           
     179 
    176180         IF ( ll_zAimp ) THEN 
    177181            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 
    … …  
    186190            DO_3D_00_00( 1, jpkm1 ) 
    187191               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
    188                   &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     192                  &                                     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    189193            END_3D 
    190194            ! 
    191195         END IF 
    192          !                 
     196         ! 
    193197         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
    194198            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:) 
    195199         END IF 
    196200         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    197          IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:)  
     201         IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:) 
    198202         ! 
    199203         !        !==  anti-diffusive flux : high order minus low order  ==! 
    … …  
    225229               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points 
    226230               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
    227                !                                                  ! C4 minus upstream advective fluxes  
     231               !                                                  ! C4 minus upstream advective fluxes 
    228232               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 
    229233               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 
    … …  
    245249               zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) ) 
    246250               zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) ) 
    247                !                                                  ! C4 minus upstream advective fluxes  
     251               !                                                  ! C4 minus upstream advective fluxes 
    248252               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 
    249253               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 
    … …  
    251255            ! 
    252256         END SELECT 
    253          !                       
     257         ! 
    254258         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values) 
    255259         ! 
    … …  
    270274            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked 
    271275         ENDIF 
    272          !          
     276         ! 
    273277         IF ( ll_zAimp ) THEN 
    274278            DO_3D_00_00( 1, jpkm1 ) 
    … …  
    277281                  &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    278282                  &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
    279                ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     283               ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs)*tmask(ji,jj,jk) 
    280284            END_3D 
    281285            ! 
    … …  
    316320            DO_3D_00_00( 1, jpkm1 ) 
    317321               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
    318                   &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    319             END_3D 
    320          END IF          
     322                  &                                     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     323            END_3D 
     324         END IF 
    321325         ! 
    322326         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport 
    323             ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes  
     327            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes 
    324328            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  !     to upstream fluxes 
    325329            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! 
    … …  
    344348         DEALLOCATE( zwdia, zwinf, zwsup ) 
    345349      ENDIF 
    346       IF( l_trd .OR. l_hst ) THEN  
     350      IF( l_trd .OR. l_hst ) THEN 
    347351         DEALLOCATE( ztrdx, ztrdy, ztrdz ) 
    348352      ENDIF 
    349       IF( l_ptr ) THEN  
     353      IF( l_ptr ) THEN 
    350354         DEALLOCATE( zptry ) 
    351355      ENDIF 
    … …  
    357361      !!--------------------------------------------------------------------- 
    358362      !!                    ***  ROUTINE nonosc  *** 
    359       !!      
    360       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
    361       !!       scheme and the before field by a nonoscillatory algorithm  
     363      !! 
     364      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
     365      !!       scheme and the before field by a nonoscillatory algorithm 
    362366      !! 
    363367      !! **  Method  :   ... ??? 
    … …  
    367371      !!       in-space based differencing for fluid 
    368372      !!---------------------------------------------------------------------- 
    369       INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index  
     373      INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index 
    370374      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step 
    371375      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field 
    … …  
    453457      !!---------------------------------------------------------------------- 
    454458      !!                  ***  ROUTINE interp_4th_cpt_org  *** 
    455       !!  
     459      !! 
    456460      !! **  Purpose :   Compute the interpolation of tracer at w-point 
    457461      !! 
    … …  
    464468      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 
    465469      !!---------------------------------------------------------------------- 
    466        
     470 
    467471      DO_3D_11_11( 3, jpkm1 ) 
    468472         zwd (ji,jj,jk) = 4._wp 
    … …  
    475479            zwi (ji,jj,jk) = 0._wp 
    476480            zws (ji,jj,jk) = 0._wp 
    477             zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )     
     481            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
    478482         ENDIF 
    479483      END_3D 
    … …  
    499503      END_2D 
    500504      DO_3D_11_11( 3, jpkm1 ) 
    501          pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
     505         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
    502506      END_3D 
    503507 
    … …  
    508512         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    509513      END_3D 
    510       !     
     514      ! 
    511515   END SUBROUTINE interp_4th_cpt_org 
    512     
     516 
    513517 
    514518   SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 
    515519      !!---------------------------------------------------------------------- 
    516520      !!                  ***  ROUTINE interp_4th_cpt  *** 
    517       !!  
     521      !! 
    518522      !! **  Purpose :   Compute the interpolation of tracer at w-point 
    519523      !! 
    … …  
    543547!      CASE( np_CEN2 )   ! 2nd order centered  at top & bottom 
    544548!      END SELECT 
    545 !!gm   
     549!!gm 
    546550      ! 
    547551      IF ( ln_isfcav ) THEN            ! set level two values which may not be set in ISF case 
    … …  
    561565         zwi (ji,jj,ikb) = 0._wp 
    562566         zws (ji,jj,ikb) = 0._wp 
    563          zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )             
     567         zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 
    564568      END_2D 
    565569      ! 
    … …  
    577581      END_2D 
    578582      DO_3D_00_00( 3, jpkm1 ) 
    579          pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
     583         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
    580584      END_3D 
    581585 
    … …  
    586590         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    587591      END_3D 
    588       !     
     592      ! 
    589593   END SUBROUTINE interp_4th_cpt 
    590594 
    … …  
    593597      !!---------------------------------------------------------------------- 
    594598      !!                  ***  ROUTINE tridia_solver  *** 
    595       !!  
     599      !! 
    596600      !! **  Purpose :   solve a symmetric 3diagonal system 
    597601      !! 
    598602      !! **  Method  :   solve M.t_out = RHS(t)  where M is a tri diagonal matrix ( jpk*jpk ) 
    599       !!      
     603      !! 
    600604      !!             ( D_1 U_1  0   0   0  )( t_1 )   ( RHS_1 ) 
    601605      !!             ( L_2 D_2 U_2  0   0  )( t_2 )   ( RHS_2 ) 
    … …  
    603607      !!             (        ...          )( ... )   ( ...  ) 
    604608      !!             (  0   0   0  L_k D_k )( t_k )   ( RHS_k ) 
    605       !!      
     609      !! 
    606610      !!        M is decomposed in the product of an upper and lower triangular matrix. 
    607       !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL  
     611      !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL 
    608612      !!        (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 
    609613      !!        The solution is pta. 
    … …  
    613617      REAL(wp),DIMENSION(:,:,:), INTENT(in   ) ::   pRHS          ! Right-Hand-Side 
    614618      REAL(wp),DIMENSION(:,:,:), INTENT(  out) ::   pt_out        !!gm field at level=F(klev) 
    615       INTEGER                  , INTENT(in   ) ::   klev          ! =1 pt_out at w-level  
     619      INTEGER                  , INTENT(in   ) ::   klev          ! =1 pt_out at w-level 
    616620      !                                                           ! =0 pt at t-level 
    617621      INTEGER ::   ji, jj, jk   ! dummy loop integers 
    … …  
    633637      END_2D 
    634638      DO_3D_00_00( kstart+1, jpkm1 ) 
    635          pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
     639         pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
    636640      END_3D 
    637641 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_mus.F90

    r12377 r12590  
    2929   USE in_out_manager ! I/O manager 
    3030   USE lib_mpp        ! distribued memory computing 
    31    USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    32    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     31   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
     32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3333 
    3434   IMPLICIT NONE 
    … …  
    3636 
    3737   PUBLIC   tra_adv_mus   ! routine called by traadv.F90 
    38     
     38 
    3939   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits 
    4040   !                                                           !  and in closed seas (orca 2 and 1 configurations) 
    4141   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index 
    42     
     42 
    4343   LOGICAL  ::   l_trd   ! flag to compute trends 
    4444   LOGICAL  ::   l_ptr   ! flag to compute poleward transport 
    … …  
    4747   !! * Substitutions 
    4848#  include "do_loop_substitute.h90" 
     49#  include "domzgr_substitute.h90" 
    4950   !!---------------------------------------------------------------------- 
    5051   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    51    !! $Id$  
     52   !! $Id$ 
    5253   !! Software governed by the CeCILL license (see ./LICENSE) 
    5354   !!---------------------------------------------------------------------- 
    … …  
    6465      !! 
    6566      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries 
    66       !!              ld_msc_ups=T :  
     67      !!              ld_msc_ups=T : 
    6768      !! 
    6869      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
    … …  
    8889      REAL(wp) ::   zv, z0v, zzwy, z0w           !   -      - 
    8990      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zslpx   ! 3D workspace 
    90       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      -  
     91      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      - 
    9192      !!---------------------------------------------------------------------- 
    9293      ! 
    … …  
    112113                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area 
    113114            END DO 
    114          ENDIF  
    115          ! 
    116       ENDIF  
    117       !       
     115         ENDIF 
     116         ! 
     117      ENDIF 
     118      ! 
    118119      l_trd = .FALSE. 
    119120      l_hst = .FALSE. 
    120121      l_ptr = .FALSE. 
    121122      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
    122       IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.  
     123      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
    123124      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    124125         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE. 
    … …  
    130131         !                                !-- first guess of the slopes 
    131132         zwx(:,:,jpk) = 0._wp                   ! bottom values 
    132          zwy(:,:,jpk) = 0._wp   
     133         zwy(:,:,jpk) = 0._wp 
    133134         DO_3D_10_10( 1, jpkm1 ) 
    134135            zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
    … …  
    176177         DO_3D_00_00( 1, jpkm1 ) 
    177178            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       & 
    178             &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
    179             &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     179            &                                           +  zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
     180            &                           * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    180181         END_3D 
    181182         !                                ! trend diagnostics 
    … …  
    184185            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) ) 
    185186         END IF 
    186          !                                 ! "Poleward" heat and salt transports  
     187         !                                 ! "Poleward" heat and salt transports 
    187188         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 
    188189         !                                 !  heat transport 
    … …  
    227228         ! 
    228229         DO_3D_00_00( 1, jpkm1 ) 
    229             pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     230            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )   & 
     231               &                                      * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    230232         END_3D 
    231233         !                                ! send trends for diagnostic 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_qck.F90

    r12377 r12590  
    1919   USE trc_oce         ! share passive tracers/Ocean variables 
    2020   USE trd_oce         ! trends: ocean variables 
    21    USE trdtra          ! trends manager: tracers  
     21   USE trdtra          ! trends manager: tracers 
    2222   USE diaptr          ! poleward transport diagnostics 
    2323   USE iom 
    … …  
    2626   USE lib_mpp         ! distribued memory computing 
    2727   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    28    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     28   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    2929 
    3030   IMPLICIT NONE 
    … …  
    4141   !! * Substitutions 
    4242#  include "do_loop_substitute.h90" 
     43#  include "domzgr_substitute.h90" 
    4344   !!---------------------------------------------------------------------- 
    4445   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    104105      l_ptr = .FALSE. 
    105106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    106       IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.  
     107      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 
    107108      ! 
    108109      ! 
    109110      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 
    110       CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )  
    111       CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )  
     111      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 
     112      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 
    112113 
    113114      !        ! vertical fluxes are computed with the 2nd order centered scheme 
    … …  
    137138      DO jn = 1, kjpt                                            ! tracer loop 
    138139         !                                                       ! =========== 
    139          zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp  
    140          zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp    
     140         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
     141         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp 
    141142         ! 
    142143!!gm why not using a SHIFT instruction... 
    … …  
    145146            zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer 
    146147         END_3D 
    147          CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions  
    148           
     148         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions 
     149 
    149150         ! 
    150151         ! Horizontal advective fluxes 
    151152         ! --------------------------- 
    152153         DO_3D_00_00( 1, jpkm1 ) 
    153             zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
    154             zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T  
    155          END_3D 
    156          ! 
    157          DO_3D_00_00( 1, jpkm1 ) 
    158             zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     154            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
     155            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T 
     156         END_3D 
     157         ! 
     158         DO_3D_00_00( 1, jpkm1 ) 
     159            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
    159160            zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
    160161            zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
    … …  
    162163            zfd(ji,jj,jk)  = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji  ,jj,jk,jn,Kbb)  ! FD in the x-direction for T 
    163164         END_3D 
    164          !--- Lateral boundary conditions  
     165         !--- Lateral boundary conditions 
    165166         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  zwx(:,:,:), 'T', 1. ) 
    166167 
    … …  
    172173            zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 
    173174         END_3D 
    174          CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions  
     175         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions 
    175176 
    176177         ! 
    177178         ! Tracer flux on the x-direction 
    178          DO jk = 1, jpkm1   
     179         DO jk = 1, jpkm1 
    179180            ! 
    180181            DO_2D_00_00 
    181                zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     182               zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
    182183               !--- If the second ustream point is a land point 
    183184               !--- the flux is computed by the 1st order UPWIND scheme 
    … …  
    226227      DO jn = 1, kjpt                                            ! tracer loop 
    227228         !                                                       ! =========== 
    228          zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0   
    229          zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0      
    230          !                                                   
    231          DO jk = 1, jpkm1                                 
    232             !                                              
     229         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
     230         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0 
     231         ! 
     232         DO jk = 1, jpkm1 
     233            ! 
    233234            !--- Computation of the ustream and downstream value of the tracer and the mask 
    234235            DO_2D_00_00 
    … …  
    239240            END_2D 
    240241         END DO 
    241          CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions  
    242  
    243           
     242         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions 
     243 
     244 
    244245         ! 
    245246         ! Horizontal advective fluxes 
    … …  
    247248         ! 
    248249         DO_3D_00_00( 1, jpkm1 ) 
    249             zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
    250             zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T  
    251          END_3D 
    252          ! 
    253          DO_3D_00_00( 1, jpkm1 ) 
    254             zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
    255             zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     250            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
     251            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T 
     252         END_3D 
     253         ! 
     254         DO_3D_00_00( 1, jpkm1 ) 
     255            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
     256            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) )   & 
     257               &         * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
    256258            zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
    257259            zfc(ji,jj,jk)  = zdir * pt(ji,jj  ,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj+1,jk,jn,Kbb)  ! FC in the x-direction for T 
    … …  
    259261         END_3D 
    260262 
    261          !--- Lateral boundary conditions  
     263         !--- Lateral boundary conditions 
    262264         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. ) 
    263265 
    … …  
    269271            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 
    270272         END_3D 
    271          CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions  
     273         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions 
    272274         ! 
    273275         ! Tracer flux on the x-direction 
    274          DO jk = 1, jpkm1   
     276         DO jk = 1, jpkm1 
    275277            ! 
    276278            DO_2D_00_00 
    277                zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     279               zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
    278280               !--- If the second ustream point is a land point 
    279281               !--- the flux is computed by the 1st order UPWIND scheme 
    … …  
    312314      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    313315      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers 
    314       REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity  
     316      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity 
    315317      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    316318      ! 
    … …  
    332334            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean) 
    333335               DO_2D_11_11 
    334                   zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface  
     336                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface 
    335337               END_2D 
    336338            ELSE                                   ! no ocean cavities (only ocean surface) 
    … …  
    356358      !! ** Purpose :  Computation of advective flux with Quickest scheme 
    357359      !! 
    358       !! ** Method :    
     360      !! ** Method : 
    359361      !!---------------------------------------------------------------------- 
    360362      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point 
    … …  
    363365      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux 
    364366      !! 
    365       INTEGER  ::  ji, jj, jk               ! dummy loop indices  
    366       REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars           
     367      INTEGER  ::  ji, jj, jk               ! dummy loop indices 
     368      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars 
    367369      REAL(wp) ::  zc, zcurv, zfho          !   -      - 
    368370      !---------------------------------------------------------------------- 
    … …  
    374376         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) ) 
    375377         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv 
    376          zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST  
     378         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST 
    377379         ! 
    378380         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk) 
    … …  
    380382         zcoef3 = ABS( zcurv ) 
    381383         IF( zcoef3 >= zcoef2 ) THEN 
    382             zfho = pfc(ji,jj,jk)  
     384            zfho = pfc(ji,jj,jk) 
    383385         ELSE 
    384386            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF 
    385387            IF( zcoef1 >= 0. ) THEN 
    386                zfho = MAX( pfc(ji,jj,jk), zfho )  
    387                zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )  
     388               zfho = MAX( pfc(ji,jj,jk), zfho ) 
     389               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
    388390            ELSE 
    389                zfho = MIN( pfc(ji,jj,jk), zfho )  
    390                zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )  
     391               zfho = MIN( pfc(ji,jj,jk), zfho ) 
     392               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
    391393            ENDIF 
    392394         ENDIF 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_ubs.F90

    r12377 r12590  
    1010   !!---------------------------------------------------------------------- 
    1111   !!   tra_adv_ubs : update the tracer trend with the horizontal 
    12    !!                 advection trends using a third order biaised scheme   
     12   !!                 advection trends using a third order biaised scheme 
    1313   !!---------------------------------------------------------------------- 
    1414   USE oce            ! ocean dynamics and active tracers 
    … …  
    1616   USE trc_oce        ! share passive tracers/Ocean variables 
    1717   USE trd_oce        ! trends: ocean variables 
    18    USE traadv_fct      ! acces to routine interp_4th_cpt  
    19    USE trdtra         ! trends manager: tracers  
     18   USE traadv_fct      ! acces to routine interp_4th_cpt 
     19   USE trdtra         ! trends manager: tracers 
    2020   USE diaptr         ! poleward transport diagnostics 
    2121   USE diaar5         ! AR5 diagnostics 
    … …  
    2525   USE lib_mpp        ! massively parallel library 
    2626   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    27    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     27   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    2828 
    2929   IMPLICIT NONE 
    … …  
    3939   !! * Substitutions 
    4040#  include "do_loop_substitute.h90" 
     41#  include "domzgr_substitute.h90" 
    4142   !!---------------------------------------------------------------------- 
    4243   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5051      !!---------------------------------------------------------------------- 
    5152      !!                  ***  ROUTINE tra_adv_ubs  *** 
    52       !!                  
     53      !! 
    5354      !! ** Purpose :   Compute the now trend due to the advection of tracers 
    5455      !!      and add it to the general trend of passive tracer equations. 
    … …  
    5960      !!      For example the i-component of the advective fluxes are given by : 
    6061      !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0 
    61       !!          ztu = !  or  
     62      !!          ztu = !  or 
    6263      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0 
    6364      !!      where zltu is the second derivative of the before temperature field: 
    6465      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 
    65       !!        This results in a dissipatively dominant (i.e. hyper-diffusive)  
    66       !!      truncation error. The overall performance of the advection scheme  
    67       !!      is similar to that reported in (Farrow and Stevens, 1995).  
     66      !!        This results in a dissipatively dominant (i.e. hyper-diffusive) 
     67      !!      truncation error. The overall performance of the advection scheme 
     68      !!      is similar to that reported in (Farrow and Stevens, 1995). 
    6869      !!        For stability reasons, the first term of the fluxes which corresponds 
    69       !!      to a second order centered scheme is evaluated using the now velocity  
    70       !!      (centered in time) while the second term which is the diffusive part  
    71       !!      of the scheme, is evaluated using the before velocity (forward in time).  
     70      !!      to a second order centered scheme is evaluated using the now velocity 
     71      !!      (centered in time) while the second term which is the diffusive part 
     72      !!      of the scheme, is evaluated using the before velocity (forward in time). 
    7273      !!      Note that UBS is not positive. Do not use it on passive tracers. 
    7374      !!                On the vertical, the advection is evaluated using a FCT scheme, 
    74       !!      as the UBS have been found to be too diffusive.  
    75       !!                kn_ubs_v argument controles whether the FCT is based on  
    76       !!      a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact  
     75      !!      as the UBS have been found to be too diffusive. 
     76      !!                kn_ubs_v argument controles whether the FCT is based on 
     77      !!      a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact 
    7778      !!      scheme (kn_ubs_v=4). 
    7879      !! 
    … …  
    8182      !!             - poleward advective heat and salt transport (ln_diaptr=T) 
    8283      !! 
    83       !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.  
    84       !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.  
     84      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 
     85      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731�1741. 
    8586      !!---------------------------------------------------------------------- 
    8687      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index 
    … …  
    111112      l_ptr = .FALSE. 
    112113      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
    113       IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.  
     114      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
    114115      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    115116         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE. 
    … …  
    122123      DO jn = 1, kjpt                                            ! tracer loop 
    123124         !                                                       ! =========== 
    124          !                                               
     125         ! 
    125126         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==! 
    126127            DO_2D_10_10 
    … …  
    135136               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
    136137            END_2D 
    137             !                                     
    138          END DO          
     138            ! 
     139         END DO 
    139140         CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1. )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    140          !     
     141         ! 
    141142         DO_3D_10_10( 1, jpkm1 ) 
    142143            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2) 
    … …  
    158159               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        & 
    159160                  &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    & 
    160                   &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     161                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) & 
     162                  &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    161163            END_2D 
    162             !                                              
     164            ! 
    163165         END DO 
    164166         ! 
    165167         zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case 
    166          !                                            ! and/or in trend diagnostic (l_trd=T)  
    167          !                 
     168         !                                            ! and/or in trend diagnostic (l_trd=T) 
     169         ! 
    168170         IF( l_trd ) THEN                  ! trend diagnostics 
    169171             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztu, pU, pt(:,:,:,jn,Kmm) ) 
    170172             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztv, pV, pt(:,:,:,jn,Kmm) ) 
    171173         END IF 
    172          !      
     174         ! 
    173175         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    174176         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) ) 
    … …  
    181183         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme 
    182184         ! 
    183          CASE(  2  )                   ! 2nd order FCT  
    184             !          
     185         CASE(  2  )                   ! 2nd order FCT 
     186            ! 
    185187            IF( l_trd )   zltv(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag. 
    186188            ! 
    … …  
    194196               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface 
    195197                  DO_2D_11_11 
    196                      ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface  
     198                     ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
    197199                  END_2D 
    198200               ELSE                                ! no cavities: only at the ocean surface 
    … …  
    202204            ! 
    203205            DO_3D_00_00( 1, jpkm1 ) 
    204                ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    205                pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak  
     206               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    & 
     207                  &     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     208               pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak 
    206209               zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
    207210            END_3D 
    … …  
    228231         ! 
    229232         DO_3D_00_00( 1, jpkm1 ) 
    230             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     233            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    & 
     234               &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    231235         END_3D 
    232236         ! 
    … …  
    235239               zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          & 
    236240                  &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   & 
    237                   &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     241                  &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    238242            END_3D 
    239243            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zltv ) 
    … …  
    248252      !!--------------------------------------------------------------------- 
    249253      !!                    ***  ROUTINE nonosc_z  *** 
    250       !!      
    251       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
    252       !!       scheme and the before field by a nonoscillatory algorithm  
     254      !! 
     255      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
     256      !!       scheme and the before field by a nonoscillatory algorithm 
    253257      !! 
    254258      !! **  Method  :   ... ??? 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trabbc.F90

    r12377 r12590  
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   tra_bbc       : update the tracer trend at ocean bottom  
     14   !!   tra_bbc       : update the tracer trend at ocean bottom 
    1515   !!   tra_bbc_init  : initialization of geothermal heat flux trend 
    1616   !!---------------------------------------------------------------------- 
    … …  
    1919   USE phycst         ! physical constants 
    2020   USE trd_oce        ! trends: ocean variables 
    21    USE trdtra         ! trends manager: tracers  
     21   USE trdtra         ! trends manager: tracers 
    2222   ! 
    2323   USE in_out_manager ! I/O manager 
    24    USE iom            ! xIOS  
     24   USE iom            ! xIOS 
    2525   USE fldread        ! read input fields 
    2626   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    … …  
    4343 
    4444   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh   ! structure of input qgh (file informations, fields read) 
    45   
     45 
    4646   !! * Substitutions 
    4747#  include "do_loop_substitute.h90" 
     48#  include "domzgr_substitute.h90" 
    4849   !!---------------------------------------------------------------------- 
    4950   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5758      !!                  ***  ROUTINE tra_bbc  *** 
    5859      !! 
    59       !! ** Purpose :   Compute the bottom boundary contition on temperature  
    60       !!              associated with geothermal heating and add it to the  
     60      !! ** Purpose :   Compute the bottom boundary contition on temperature 
     61      !!              associated with geothermal heating and add it to the 
    6162      !!              general trend of temperature equations. 
    6263      !! 
    63       !! ** Method  :   The geothermal heat flux set to its constant value of  
     64      !! ** Method  :   The geothermal heat flux set to its constant value of 
    6465      !!              86.4 mW/m2 (Stein and Stein 1992, Huang 1999). 
    6566      !!       The temperature trend associated to this heat flux through the 
    … …  
    9192      !                             !  Add the geothermal trend on temperature 
    9293      DO_2D_00_00 
    93          pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm) 
     94         pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs)   & 
     95            &             + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm) 
    9496      END_2D 
    9597      ! 
    … …  
    133135      CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files 
    134136      !! 
    135       NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir  
     137      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 
    136138      !!---------------------------------------------------------------------- 
    137139      ! 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trabbl.F90

    r12377 r12590  
    3131   USE trdtra         ! trends: active tracers 
    3232   ! 
    33    USE iom            ! IOM library                
     33   USE iom            ! IOM library 
    3434   USE in_out_manager ! I/O manager 
    3535   USE lbclnk         ! ocean lateral boundary conditions 
    3636   USE prtctl         ! Print control 
    3737   USE timing         ! Timing 
    38    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     38   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3939 
    4040   IMPLICIT NONE 
    … …  
    6868   !! * Substitutions 
    6969#  include "do_loop_substitute.h90" 
     70#  include "domzgr_substitute.h90" 
    7071   !!---------------------------------------------------------------------- 
    7172   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    195196            zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S 
    196197         END_2D 
    197          !                
     198         ! 
    198199         DO_2D_00_00 
    199200            ik = mbkt(ji,jj)                            ! bottom T-level index 
    … …  
    391392               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
    392393               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
    393                !                                                          ! 2*masked bottom density gradient  
     394               !                                                          ! 2*masked bottom density gradient 
    394395               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
    395396                         - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
    … …  
    513514      END_2D 
    514515      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
    515       zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp )   
    516       CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1., zmbkv,'V',1.)  
     516      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) 
     517      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1., zmbkv,'V',1.) 
    517518      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 ) 
    518519      ! 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traisf.F90

    r12377 r12590  
    3333      !!---------------------------------------------------------------------- 
    3434      !!                  ***  ROUTINE tra_isf  *** 
    35       !!                    
     35      !! 
    3636      !! ** Purpose :  Compute the temperature trend due to the ice shelf melting (qhoce + qhc) 
    3737      !! 
    … …  
    6161         ! 
    6262         ! Dynamical stability at start up after change in under ice shelf cavity geometry is achieve by correcting the divergence. 
    63          ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping  
     63         ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping 
    6464         ! the same as at the end of the latest time step. So correction need to be apply at nit000 (euler time step) and 
    6565         ! half of it at nit000+1 (leap frog time step). 
    … …  
    8989      !! *** Purpose :  Compute the temperature trend due to the ice shelf melting (qhoce + qhc) for cav or par case 
    9090      !! 
    91       !! *** Action :: Update pts(:,:,:,:,Krhs) with the surface boundary condition trend  
     91      !! *** Action :: Update pts(:,:,:,:,Krhs) with the surface boundary condition trend 
    9292      !! 
    9393      !!---------------------------------------------------------------------- 
    … …  
    9898      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc , ptsc_b 
    9999      !!---------------------------------------------------------------------- 
    100       INTEGER                      :: ji,jj,jk  ! loop index    
     100      INTEGER                      :: ji,jj,jk  ! loop index 
    101101      INTEGER                      :: ikt, ikb  ! top and bottom level of the tbl 
    102102      REAL(wp), DIMENSION(jpi,jpj) :: ztc       ! total ice shelf tracer trend 
    … …  
    117117         END DO 
    118118         ! 
    119          ! level partially include in ice shelf boundary layer  
     119         ! level partially include in ice shelf boundary layer 
    120120         pts(ji,jj,ikb,jp_tem) = pts(ji,jj,ikb,jp_tem) + ztc(ji,jj) * pfrac(ji,jj) 
    121121         ! 
    … …  
    128128      !!                  ***  ROUTINE tra_isf_cpl  *** 
    129129      !! 
    130       !! *** Action :: Update pts(:,:,:,:,Krhs) with the ice shelf coupling trend  
     130      !! *** Action :: Update pts(:,:,:,:,Krhs) with the ice shelf coupling trend 
    131131      !! 
    132132      !!---------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traldf_iso.F90

    r12377 r12590  
    1515   !!---------------------------------------------------------------------- 
    1616   !!   tra_ldf_iso   : update the tracer trend with the horizontal component of a iso-neutral laplacian operator 
    17    !!                   and with the vertical part of the isopycnal or geopotential s-coord. operator  
     17   !!                   and with the vertical part of the isopycnal or geopotential s-coord. operator 
    1818   !!---------------------------------------------------------------------- 
    1919   USE oce            ! ocean dynamics and active tracers 
    … …  
    4141   !! * Substitutions 
    4242#  include "do_loop_substitute.h90" 
     43#  include "domzgr_substitute.h90" 
    4344   !!---------------------------------------------------------------------- 
    4445   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5455      !!                  ***  ROUTINE tra_ldf_iso  *** 
    5556      !! 
    56       !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
    57       !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and  
     57      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
     58      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 
    5859      !!      add it to the general trend of tracer equation. 
    5960      !! 
    60       !! ** Method  :   The horizontal component of the lateral diffusive trends  
     61      !! ** Method  :   The horizontal component of the lateral diffusive trends 
    6162      !!      is provided by a 2nd order operator rotated along neural or geopo- 
    6263      !!      tential surfaces to which an eddy induced advection can be added 
    … …  
    6970      !! 
    7071      !!      2nd part :  horizontal fluxes of the lateral mixing operator 
    71       !!      ========     
     72      !!      ======== 
    7273      !!         zftu =  pahu e2u*e3u/e1u di[ tb ] 
    7374      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ] 
    … …  
    111112      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      - 
    112113      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d 
    113       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw  
     114      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw 
    114115      !!---------------------------------------------------------------------- 
    115116      ! 
    … …  
    119120         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    120121         ! 
    121          akz     (:,:,:) = 0._wp       
     122         akz     (:,:,:) = 0._wp 
    122123         ah_wslp2(:,:,:) = 0._wp 
    123124      ENDIF 
    124       !    
     125      ! 
    125126      l_hst = .FALSE. 
    126127      l_ptr = .FALSE. 
    127       IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE.  
     128      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE. 
    128129      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    129130         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
    … …  
    138139      ELSE                    ;   zsign = -1._wp 
    139140      ENDIF 
    140           
     141 
    141142      !!---------------------------------------------------------------------- 
    142143      !!   0 - calculate  ah_wslp2 and akz 
    … …  
    173174               DO_3D_10_10( 2, jpkm1 ) 
    174175                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    175                      &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  ) 
     176                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk)    & 
     177                     &                        / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  ) 
    176178               END_3D 
    177179            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
    … …  
    184186           ! 
    185187         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
    186             akz(:,:,:) = ah_wslp2(:,:,:)       
     188            akz(:,:,:) = ah_wslp2(:,:,:) 
    187189         ENDIF 
    188190      ENDIF 
    … …  
    191193      DO jn = 1, kjpt                                            ! tracer loop 
    192194         !                                                       ! =========== 
    193          !                                                
    194          !!---------------------------------------------------------------------- 
    195          !!   I - masked horizontal derivative  
     195         ! 
     196         !!---------------------------------------------------------------------- 
     197         !!   I - masked horizontal derivative 
    196198         !!---------------------------------------------------------------------- 
    197199!!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient.... 
    … …  
    200202         !!end 
    201203 
    202          ! Horizontal tracer gradient  
     204         ! Horizontal tracer gradient 
    203205         DO_3D_10_10( 1, jpkm1 ) 
    204206            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    … …  
    207209         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient 
    208210            DO_2D_10_10 
    209                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     211               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
    210212               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    211213            END_2D 
    212214            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity 
    213215               DO_2D_10_10 
    214                   IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
    215                   IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
     216                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 
     217                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 
    216218               END_2D 
    217219            ENDIF 
    … …  
    248250               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    249251                  &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
    250                   &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                   
     252                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
    251253            END_2D 
    252254            ! 
    … …  
    254256               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
    255257                  &                                                 + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
    256                   &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     258                  &                                             * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    257259            END_2D 
    258          END DO                                        !   End of slab   
     260         END DO                                        !   End of slab 
    259261 
    260262         !!---------------------------------------------------------------------- 
    … …  
    266268         !                          ! Surface and bottom vertical fluxes set to zero 
    267269         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
    268           
     270 
    269271         DO_3D_00_00( 2, jpkm1 ) 
    270272            ! 
    … …  
    295297            END_3D 
    296298            ! 
    297          ELSE                                   ! bilaplacian  
     299         ELSE                                   ! bilaplacian 
    298300            SELECT CASE( kpass ) 
    299301            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    … …  
    311313            END SELECT 
    312314         ENDIF 
    313          !          
     315         ! 
    314316         DO_3D_00_00( 1, jpkm1 ) 
    315317            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traldf_lap_blp.F90

    r12377 r12590  
    44   !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian) 
    55   !!============================================================================== 
    6    !! History :  3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian  
     6   !! History :  3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian 
    77   !!---------------------------------------------------------------------- 
    88 
    … …  
    3838   !! * Substitutions 
    3939#  include "do_loop_substitute.h90" 
     40#  include "domzgr_substitute.h90" 
    4041   !!---------------------------------------------------------------------- 
    4142   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    4748   SUBROUTINE tra_ldf_lap( kt, Kmm, kit000, cdtype, pahu, pahv  ,               & 
    4849      &                                             pgu , pgv   , pgui, pgvi,   & 
    49       &                                             pt  , pt_rhs, kjpt, kpass )  
     50      &                                             pt  , pt_rhs, kjpt, kpass ) 
    5051      !!---------------------------------------------------------------------- 
    5152      !!                  ***  ROUTINE tra_ldf_lap  *** 
    52       !!                    
    53       !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
     53      !! 
     54      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
    5455      !!      trend and add it to the general trend of tracer equation. 
    5556      !! 
    5657      !! ** Method  :   Second order diffusive operator evaluated using before 
    57       !!      fields (forward time scheme). The horizontal diffusive trends of  
     58      !!      fields (forward time scheme). The horizontal diffusive trends of 
    5859      !!      the tracer is given by: 
    5960      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ] 
    … …  
    6263      !!          pt_rhs = pt_rhs + difft 
    6364      !! 
    64       !! ** Action  : - Update pt_rhs arrays with the before iso-level  
     65      !! ** Action  : - Update pt_rhs arrays with the before iso-level 
    6566      !!                harmonic mixing trend. 
    6667      !!---------------------------------------------------------------------- 
    … …  
    7576      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    7677      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! before tracer fields 
    77       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend  
     78      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend 
    7879      ! 
    7980      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    105106      !                             ! =========== ! 
    106107      DO jn = 1, kjpt               ! tracer loop ! 
    107          !                          ! =========== !     
    108          !                                
     108         !                          ! =========== ! 
     109         ! 
    109110         DO_3D_10_10( 1, jpkm1 ) 
    110111            ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) 
    … …  
    118119            IF( ln_isfcav ) THEN                ! top in ocean cavities only 
    119120               DO_2D_10_10 
    120                   IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn)  
    121                   IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn)  
     121                  IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
     122                  IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 
    122123               END_2D 
    123124            ENDIF 
    … …  
    142143      ! 
    143144   END SUBROUTINE tra_ldf_lap 
    144     
     145 
    145146 
    146147   SUBROUTINE tra_ldf_blp( kt, Kmm, kit000, cdtype, pahu, pahv  ,             & 
    … …  
    149150      !!---------------------------------------------------------------------- 
    150151      !!                 ***  ROUTINE tra_ldf_blp  *** 
    151       !!                     
    152       !! ** Purpose :   Compute the before lateral tracer diffusive  
     152      !! 
     153      !! ** Purpose :   Compute the before lateral tracer diffusive 
    153154      !!      trend and add it to the general trend of tracer equation. 
    154155      !! 
    … …  
    200201      ! 
    201202      CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1. )     ! Lateral boundary conditions (unchanged sign) 
    202       !                                               ! Partial top/bottom cell: GRADh( zlap )   
     203      !                                               ! Partial top/bottom cell: GRADh( zlap ) 
    203204      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom 
    204       ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, Kmm, kjpt, zlap, zglu, zglv )              ! only bottom  
     205      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, Kmm, kjpt, zlap, zglu, zglv )              ! only bottom 
    205206      ENDIF 
    206207      ! 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traldf_triad.F90

    r12377 r12590  
    4141   !! * Substitutions 
    4242#  include "do_loop_substitute.h90" 
     43#  include "domzgr_substitute.h90" 
    4344   !!---------------------------------------------------------------------- 
    4445   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    108109         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 
    109110      ENDIF 
    110       !    
     111      ! 
    111112      l_hst = .FALSE. 
    112113      l_ptr = .FALSE. 
    113       IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )      l_ptr = .TRUE.  
     114      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )      l_ptr = .TRUE. 
    114115      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    115116         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
    … …  
    124125      ELSE                    ;   zsign = -1._wp 
    125126      ENDIF 
    126       !     
     127      ! 
    127128      !!---------------------------------------------------------------------- 
    128129      !!   0 - calculate  ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw 
    … …  
    131132      IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==! 
    132133         ! 
    133          akz     (:,:,:) = 0._wp       
     134         akz     (:,:,:) = 0._wp 
    134135         ah_wslp2(:,:,:) = 0._wp 
    135136         IF( ln_ldfeiv_dia ) THEN 
    … …  
    158159         END DO 
    159160         ! 
    160          DO jp = 0, 1                            ! j-k triads  
     161         DO jp = 0, 1                            ! j-k triads 
    161162            DO kp = 0, 1 
    162163               DO_3D_10_10( 1, jpkm1 ) 
    … …  
    184185               DO_3D_10_10( 2, jpkm1 ) 
    185186                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    186                      &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  ) 
     187                     &          * (  akz(ji,jj,jk)              & 
     188                     &            + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  ) 
    187189               END_3D 
    188190            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
    … …  
    195197           ! 
    196198         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
    197             akz(:,:,:) = ah_wslp2(:,:,:)       
     199            akz(:,:,:) = ah_wslp2(:,:,:) 
    198200         ENDIF 
    199201         ! 
    … …  
    222224            IF( ln_isfcav ) THEN                   ! top level (ocean cavities only) 
    223225               DO_2D_10_10 
    224                   IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn)  
    225                   IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn)  
     226                  IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 
     227                  IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 
    226228               END_2D 
    227229            ENDIF 
    … …  
    330332            !                             !==  horizontal divergence and add to the general trend  ==! 
    331333            DO_2D_00_00 
    332                pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       & 
    333                   &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   & 
    334                   &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  ) 
     334               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    & 
     335                  &                            + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       & 
     336                  &                                       + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   & 
     337                  &                                     / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  ) 
    335338            END_2D 
    336339            ! 
    … …  
    344347                  &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    345348            END_3D 
    346          ELSE                                   ! bilaplacian  
     349         ELSE                                   ! bilaplacian 
    347350            SELECT CASE( kpass ) 
    348351            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    … …  
    357360                     &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   ) 
    358361               END_3D 
    359             END SELECT  
     362            END SELECT 
    360363         ENDIF 
    361364         ! 
    362365         DO_3D_00_00( 1, jpkm1 ) 
    363             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
    364                &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 
     366            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    & 
     367            &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
     368               &                                       / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 
    365369         END_3D 
    366370         ! 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/tramle.F90

    r12377 r12590  
    4949   !! * Substitutions 
    5050#  include "do_loop_substitute.h90" 
     51#  include "domzgr_substitute.h90" 
    5152   !!---------------------------------------------------------------------- 
    5253   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/tranpc.F90

    r12377 r12590  
    3535   !! * Substitutions 
    3636#  include "do_loop_substitute.h90" 
     37#  include "domzgr_substitute.h90" 
    3738   !!---------------------------------------------------------------------- 
    3839   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    7172      REAL(wp), DIMENSION(        jpk     )   ::   zvn2         ! vertical profile of N2 at 1 given point... 
    7273      REAL(wp), DIMENSION(        jpk,jpts)   ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point 
    73       REAL(wp), DIMENSION(jpi,jpj,jpk     )   ::   zn2          ! N^2  
     74      REAL(wp), DIMENSION(jpi,jpj,jpk     )   ::   zn2          ! N^2 
    7475      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)   ::   zab          ! alpha and beta 
    7576      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds ! 3D workspace 
    … …  
    8687         IF( l_trdtra )   THEN                    !* Save initial after fields 
    8788            ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    88             ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa)  
     89            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 
    8990            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 
    9091         ENDIF 
    … …  
    9293         IF( l_LB_debug ) THEN 
    9394            ! Location of 1 known convection site to follow what's happening in the water column 
    94             ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...            
     95            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column... 
    9596            nncpu = 1  ;            ! the CPU domain contains the convection spot 
    96             klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...           
     97            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point... 
    9798         ENDIF 
    9899         ! 
    … …  
    105106            ! 
    106107            IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points 
    107                !                                     ! consider one ocean column  
     108               !                                     ! consider one ocean column 
    108109               zvts(:,jp_tem) = pts(ji,jj,:,jp_tem,Kaa)      ! temperature 
    109110               zvts(:,jp_sal) = pts(ji,jj,:,jp_sal,Kaa)      ! salinity 
    110111               ! 
    111                zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha  
    112                zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta   
    113                zvn2(:)         = zn2(ji,jj,:)            ! N^2  
     112               zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha 
     113               zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta 
     114               zvn2(:)         = zn2(ji,jj,:)            ! N^2 
    114115               ! 
    115116               IF( l_LB_debug ) THEN                  !LB debug: 
    … …  
    117118                  IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 
    118119                  ! writing only if on CPU domain where conv region is: 
    119                   lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                       
     120                  lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 
    120121               ENDIF                                  !LB debug  end 
    121122               ! 
    … …  
    129130                  ! 
    130131                  jiter = jiter + 1 
    131                   !  
     132                  ! 
    132133                  IF( jiter >= 400 ) EXIT 
    133134                  ! 
    … …  
    144145                        ilayer = ilayer + 1    ! yet another instable portion of the water column found.... 
    145146                        ! 
    146                         IF( lp_monitor_point ) THEN  
     147                        IF( lp_monitor_point ) THEN 
    147148                           WRITE(numout,*) 
    148149                           IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability 
    … …  
    159160                        ENDIF 
    160161                        ! 
    161                         IF( jiter == 1 )   inpcc = inpcc + 1  
     162                        IF( jiter == 1 )   inpcc = inpcc + 1 
    162163                        ! 
    163164                        IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer 
    … …  
    184185                        zsum_beta = 0._wp 
    185186                        zsum_z    = 0._wp 
    186                                                   
     187 
    187188                        DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column 
    188189                           ! 
    … …  
    193194                           zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 
    194195                           zsum_z    = zsum_z    + zdz 
    195                            !                               
     196                           ! 
    196197                           IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line 
    197198                           !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0): 
    198199                           IF( zvn2(jk+1) > zn2_zero ) EXIT 
    199200                        END DO 
    200                         
     201 
    201202                        ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2 
    202203                        IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2') 
    … …  
    224225                           zvab(jk,jp_sal) = zbeta 
    225226                        END DO 
    226                          
    227                          
     227 
     228 
    228229                        !! Updating N2 in the relvant portion of the water column 
    229230                        !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
    230231                        !! => Need to re-compute N2! will use Alpha and Beta! 
    231                          
     232 
    232233                        ikup   = MAX(2,ikup)         ! ikup can never be 1 ! 
    233234                        ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown! 
    234                          
     235 
    235236                        DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown! 
    236237 
    … …  
    252253 
    253254                        END DO 
    254                       
     255 
    255256                        ikp = MIN(ikdown+1,ikbot) 
    256                          
     257 
    257258 
    258259                     ENDIF  !IF( zvn2(ikp) < 0. ) 
    … …  
    264265 
    265266                  IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3') 
    266                   
     267 
    267268                  ! ******* At this stage ikp == ikbot ! ******* 
    268                   
     269 
    269270                  IF( ilayer > 0 ) THEN      !! least an unstable layer has been found 
    270271                     ! 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traqsr.F90

    r12377 r12590  
    99   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
    1010   !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate 
    11    !!            3.2  !  2009-04  (G. Madec & NEMO team)  
    12    !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model  
     11   !!            3.2  !  2009-04  (G. Madec & NEMO team) 
     12   !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
    1313   !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 
    14    !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume  
     14   !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume 
    1515   !!---------------------------------------------------------------------- 
    1616 
    1717   !!---------------------------------------------------------------------- 
    18    !!   tra_qsr       : temperature trend due to the penetration of solar radiation  
    19    !!   tra_qsr_init  : initialization of the qsr penetration  
     18   !!   tra_qsr       : temperature trend due to the penetration of solar radiation 
     19   !!   tra_qsr_init  : initialization of the qsr penetration 
    2020   !!---------------------------------------------------------------------- 
    2121   USE oce            ! ocean dynamics and active tracers 
    … …  
    4444   !                                 !!* Namelist namtra_qsr: penetrative solar radiation 
    4545   LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag 
    46    LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag   
     46   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag 
    4747   LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag 
    4848   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag 
    … …  
    5353   ! 
    5454   INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m) 
    55   
     55 
    5656   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll 
    5757   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data 
    … …  
    6868   !! * Substitutions 
    6969#  include "do_loop_substitute.h90" 
     70#  include "domzgr_substitute.h90" 
    7071   !!---------------------------------------------------------------------- 
    7172   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    8687      !!      Considering the 2 wavebands case: 
    8788      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 
    88       !!         The temperature trend associated with the solar radiation penetration  
     89      !!         The temperature trend associated with the solar radiation penetration 
    8990      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp) 
    9091      !!         At the bottom, boudary condition for the radiation is no flux : 
    9192      !!      all heat which has not been absorbed in the above levels is put 
    9293      !!      in the last ocean level. 
    93       !!         The computation is only done down to the level where  
    94       !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) .  
     94      !!         The computation is only done down to the level where 
     95      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) . 
    9596      !! 
    9697      !! ** Action  : - update ta with the penetrative solar radiation trend 
    … …  
    112113      REAL(wp) ::   zz0 , zz1                !    -         - 
    113114      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 
    114       REAL(wp) ::   zlogc, zlogc2, zlogc3  
     115      REAL(wp) ::   zlogc, zlogc2, zlogc3 
    115116      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr 
    116117      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 
    … …  
    127128      ! 
    128129      IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend 
    129          ALLOCATE( ztrdt(jpi,jpj,jpk) )  
     130         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
    130131         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 
    131132      ENDIF 
    … …  
    163164         ALLOCATE( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , & 
    164165            &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , & 
    165             &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   )  
     166            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   ) 
    166167         ! 
    167168         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll 
    … …  
    183184                     zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 
    184185                     zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 
    185                      zCze    = 1.12  * (zchl)**0.803  
     186                     zCze    = 1.12  * (zchl)**0.803 
    186187                     ! 
    187188                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) ) 
    … …  
    192193         ELSE                                !* constant chrlorophyll 
    193194           DO jk = 1, nksr + 1 
    194               zchl3d(:,:,jk) = 0.05  
     195              zchl3d(:,:,jk) = 0.05 
    195196            ENDDO 
    196197         ENDIF 
    … …  
    231232         END_3D 
    232233         ! 
    233          DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d )  
     234         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
    234235         ! 
    235236      CASE( np_2BD  )            !==  2-bands fluxes  ==! 
    … …  
    240241            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r ) 
    241242            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) 
    242             qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) )  
     243            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
    243244         END_3D 
    244245         ! 
    … …  
    248249      DO_3D_00_00( 1, nksr ) 
    249250         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
    250             &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm) 
     251            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) )   & 
     252            &                             / e3t(ji,jj,jk,Kmm) 
    251253      END_3D 
    252254      ! 
    … …  
    264266         DO jk = nksr, 1, -1 
    265267            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rau0_rcp 
    266          END DO          
     268         END DO 
    267269         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation 
    268          DEALLOCATE( zetot )  
     270         DEALLOCATE( zetot ) 
    269271      ENDIF 
    270272      ! 
    … …  
    272274         IF( lwxios ) CALL iom_swap(      cwxios_context          ) 
    273275         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios ) 
    274          CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios )  
     276         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
    275277         IF( lwxios ) CALL iom_swap(      cxios_context          ) 
    276278      ENDIF 
    … …  
    279281         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
    280282         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
    281          DEALLOCATE( ztrdt )  
     283         DEALLOCATE( ztrdt ) 
    282284      ENDIF 
    283285      !                       ! print mean trends (used for debugging) 
    … …  
    298300      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio 
    299301      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The 
    300       !!      default values correspond to clear water (type I in Jerlov'  
     302      !!      default values correspond to clear water (type I in Jerlov' 
    301303      !!      (1968) classification. 
    302304      !!         called by tra_qsr at the first timestep (nit000) 
    … …  
    348350         &                               ' 2 bands, 3 RGB bands or bio-model light penetration' ) 
    349351      ! 
    350       IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB  
     352      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
    351353      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc 
    352354      IF( ln_qsr_2bd                      )   nqsr = np_2BD 
    … …  
    358360      ! 
    359361      SELECT CASE( nqsr ) 
    360       !                                
     362      ! 
    361363      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==! 
    362          !                              
     364         ! 
    363365         IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration ' 
    364366         ! 
    365367         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef. 
    366          !                                    
     368         ! 
    367369         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction 
    368370         ! 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trasbc.F90

    r12377 r12590  
    99   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    1010   !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
    11    !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing  
     11   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing 
    1212   !!            4.1  !  2019-09  (P. Mathiot) isf moved in traisf 
    1313   !!---------------------------------------------------------------------- 
    … …  
    2121   USE phycst         ! physical constant 
    2222   USE eosbn2         ! Equation Of State 
    23    USE sbcmod         ! ln_rnf   
    24    USE sbcrnf         ! River runoff   
     23   USE sbcmod         ! ln_rnf 
     24   USE sbcrnf         ! River runoff 
    2525   USE traqsr         ! solar radiation penetration 
    2626   USE trd_oce        ! trends: ocean variables 
    27    USE trdtra         ! trends manager: tracers  
    28 #if defined key_asminc    
     27   USE trdtra         ! trends manager: tracers 
     28#if defined key_asminc 
    2929   USE asminc         ! Assimilation increment 
    3030#endif 
    … …  
    4343   !! * Substitutions 
    4444#  include "do_loop_substitute.h90" 
     45#  include "domzgr_substitute.h90" 
    4546   !!---------------------------------------------------------------------- 
    4647   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5354      !!---------------------------------------------------------------------- 
    5455      !!                  ***  ROUTINE tra_sbc  *** 
    55       !!                    
     56      !! 
    5657      !! ** Purpose :   Compute the tracer surface boundary condition trend of 
    5758      !!      (flux through the interface, concentration/dilution effect) 
    5859      !!      and add it to the general trend of tracer equations. 
    5960      !! 
    60       !! ** Method :   The (air+ice)-sea flux has two components:  
    61       !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface);  
    62       !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice.  
     61      !! ** Method :   The (air+ice)-sea flux has two components: 
     62      !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface); 
     63      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice. 
    6364      !!               The input forcing fields (emp, rnf, sfx) contain Fext+Fwe, 
    6465      !!             they are simply added to the tracer trend (ts(Krhs)). 
    … …  
    6869      !!             concentration/dilution effect associated with water exchanges. 
    6970      !! 
    70       !! ** Action  : - Update ts(Krhs) with the surface boundary condition trend  
     71      !! ** Action  : - Update ts(Krhs) with the surface boundary condition trend 
    7172      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T) 
    7273      !!---------------------------------------------------------------------- 
    … …  
    7576      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts        ! active tracers and RHS of tracer equation 
    7677      ! 
    77       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices   
     78      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
    7879      INTEGER  ::   ikt, ikb                    ! local integers 
    7980      REAL(wp) ::   zfact, z1_e3t, zdep, ztim   ! local scalar 
    … …  
    9091      ! 
    9192      IF( l_trdtra ) THEN                    !* Save ta and sa trends 
    92          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )  
     93         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    9394         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 
    9495         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 
    … …  
    127128         sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting 
    128129      END_2D 
    129       IF( ln_linssh ) THEN                !* linear free surface   
     130      IF( ln_linssh ) THEN                !* linear free surface 
    130131         DO_2D_01_00 
    131132            sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) 
    … …  
    138139      DO jn = 1, jpts               !==  update tracer trend  ==! 
    139140         DO_2D_01_00 
    140             pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,Kmm) 
     141            pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) )    & 
     142               &                                                / e3t(ji,jj,1,Kmm) 
    141143         END_2D 
    142144      END DO 
    143       !                   
     145      ! 
    144146      IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==! 
    145147         IF( lwxios ) CALL iom_swap(      cwxios_context          ) 
    … …  
    153155      !---------------------------------------- 
    154156      ! 
    155       IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff  
     157      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff 
    156158         zfact = 0.5_wp 
    157159         DO_2D_01_00 
    … …  
    162164                                        &                      +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 
    163165                  IF( ln_rnf_sal )   pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)                                  & 
    164                                         &                      +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep  
     166                                        &                      +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
    165167               END DO 
    166168            ENDIF 
    … …  
    179181      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation 
    180182          ! 
    181          IF( ln_linssh ) THEN  
     183         IF( ln_linssh ) THEN 
    182184            DO_2D_01_00 
    183185               ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm) 
    … …  
    202204         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt ) 
    203205         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds ) 
    204          DEALLOCATE( ztrdt , ztrds )  
     206         DEALLOCATE( ztrdt , ztrds ) 
    205207      ENDIF 
    206208      ! 

  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trazdf.F90

    r12377 r12590  
    1313   !!---------------------------------------------------------------------- 
    1414   USE oce            ! ocean dynamics and tracers variables 
    15    USE dom_oce        ! ocean space and time domain variables  
     15   USE dom_oce        ! ocean space and time domain variables 
    1616   USE domvvl         ! variable volume 
    1717   USE phycst         ! physical constant 
    … …  
    1919   USE sbc_oce        ! surface boundary condition: ocean 
    2020   USE ldftra         ! lateral diffusion: eddy diffusivity 
    21    USE ldfslp         ! lateral diffusion: iso-neutral slope  
     21   USE ldfslp         ! lateral diffusion: iso-neutral slope 
    2222   USE trd_oce        ! trends: ocean variables 
    2323   USE trdtra         ! trends: tracer trend manager 
    … …  
    3737   !! * Substitutions 
    3838#  include "do_loop_substitute.h90" 
     39#  include "domzgr_substitute.h90" 
    3940   !!---------------------------------------------------------------------- 
    4041   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    7778      ! 
    7879      !                                      !* compute lateral mixing trend and add it to the general trend 
    79       CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts )  
     80      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 
    8081 
    8182!!gm WHY here !   and I don't like that ! 
    … …  
    8889      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    8990         DO jk = 1, jpkm1 
    90             ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 
    91                &          / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 
    92             ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 
    93               &           / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 
     91            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa)    & 
     92               &              - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) )  & 
     93               &             / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 
     94            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa)    & 
     95               &              - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 
     96               &             / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 
    9497         END DO 
    9598!!gm this should be moved in trdtra.F90 and done on all trends 
    … …  
    108111   END SUBROUTINE tra_zdf 
    109112 
    110   
    111    SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt )  
     113 
     114   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 
    112115      !!---------------------------------------------------------------------- 
    113116      !!                  ***  ROUTINE tra_zdf_imp  *** 
    114117      !! 
    115118      !! ** Purpose :   Compute the after tracer through a implicit computation 
    116       !!     of the vertical tracer diffusion (including the vertical component  
    117       !!     of lateral mixing (only for 2nd order operator, for fourth order  
    118       !!     it is already computed and add to the general trend in traldf)  
     119      !!     of the vertical tracer diffusion (including the vertical component 
     120      !!     of lateral mixing (only for 2nd order operator, for fourth order 
     121      !!     it is already computed and add to the general trend in traldf) 
    119122      !! 
    120123      !! ** Method  :  The vertical diffusion of a tracer ,t , is given by: 
    … …  
    158161            zwt(:,:,1) = 0._wp 
    159162            ! 
    160             IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution  
    161                IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator  
     163            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution 
     164               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator 
    162165                  DO_3D_00_00( 2, jpkm1 ) 
    163                      zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
     166                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 
    164167                  END_3D 
    165168               ELSE                          ! standard or triad iso-neutral operator 
    … …  
    204207            !   The solution will be in the 4d array pta. 
    205208            !   The 3d array zwt is used as a work space array. 
    206             !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then  
     209            !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 
    207210            !   used as a work space array: its value is modified. 
    208211            ! 
    … …  
    214217            END_3D 
    215218            ! 
    216          ENDIF  
    217          !          
     219         ENDIF 
     220         ! 
    218221         DO_2D_00_00 
    219             pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 
     222            pt(ji,jj,1,jn,Kaa) =        e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb)    & 
     223               &               + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 
    220224         END_2D 
    221225         DO_3D_00_00( 2, jpkm1 ) 
    222             zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side 
     226            zrhs =        e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb)    &  
     227               & + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side 
    223228            pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 
    224229         END_3D 
Note: See TracChangeset for help on using the changeset viewer.