Changeset 5770


Ignore:
Timestamp:
2015-10-01T14:48:08+02:00 (5 years ago)
Author:
gm
Message:

#1593: LDF-ADV, step II.2: phasing the improvements/simplifications of advective tracer trend (see wiki page)

Location:
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO
Files:
1 added
3 deleted
25 edited
2 moved

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5753 r5770  
    2121 
    2222   IMPLICIT NONE 
    23    PUBLIC             ! allows the acces to par_oce when dom_oce is used 
    24    !                  ! exception to coding rules... to be suppressed ??? 
     23   PUBLIC             ! allows the acces to par_oce when dom_oce is used (exception to coding rules) 
    2524 
    2625   PUBLIC dom_oce_alloc  ! Called from nemogcm.F90 
     
    107106   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdttra  !: vertical profile of tracer time step 
    108107   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dtra  !: = 2*rdttra except at nit000 (=rdttra) if neuler=0 
    109  
    110    !                                         !!* Namelist namcla : cross land advection 
    111    INTEGER, PUBLIC ::   nn_cla               !: =1 cross land advection for exchanges through some straits (ORCA2) 
    112108 
    113109   !!---------------------------------------------------------------------- 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5363 r5770  
    8181      ENDIF 
    8282      ! 
    83                              CALL dom_nam      ! read namelist ( namrun, namdom, namcla ) 
     83                             CALL dom_nam      ! read namelist ( namrun, namdom ) 
    8484                             CALL dom_clo      ! Closed seas and lake 
    8585                             CALL dom_hgr      ! Horizontal mesh 
     
    131131      !! ** input   : - namrun namelist 
    132132      !!              - namdom namelist 
    133       !!              - namcla namelist 
    134133      !!              - namnc4 namelist   ! "key_netcdf4" only 
    135134      !!---------------------------------------------------------------------- 
     
    146145         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 
    147146         &             ppa2, ppkth2, ppacr2 
    148       NAMELIST/namcla/ nn_cla 
    149147#if defined key_netcdf4 
    150148      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
     
    305303      rdth      = rn_rdth 
    306304 
    307       REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection 
    308       READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905) 
    309 905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp ) 
    310  
    311       REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection 
    312       READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 ) 
    313 906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp ) 
    314       IF(lwm) WRITE( numond, namcla ) 
    315  
    316       IF(lwp) THEN 
    317          WRITE(numout,*) 
    318          WRITE(numout,*) '   Namelist namcla' 
    319          WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla 
    320       ENDIF 
    321       IF ( nn_cla .EQ. 1 ) THEN 
    322          IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2  
    323             CONTINUE 
    324          ELSE 
    325             CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' ) 
    326          ENDIF 
    327       ENDIF 
    328  
    329305#if defined key_netcdf4 
    330306      !                             ! NetCDF 4 case   ("key_netcdf4" defined) 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5552 r5770  
    1717   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization 
    1818   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option 
     19   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    129130      !!               tmask_i  : interior ocean mask 
    130131      !!---------------------------------------------------------------------- 
    131       ! 
    132       INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     132      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
    133133      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers 
    134134      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       - 
     
    284284      ! 3. Ocean/land mask at wu-, wv- and w points  
    285285      !---------------------------------------------- 
    286       wmask (:,:,1) = tmask(:,:,1) ! ???????? 
    287       wumask(:,:,1) = umask(:,:,1) ! ???????? 
    288       wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
    289       DO jk=2,jpk 
    290          wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 
    291          wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)    
    292          wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 
     286      wmask (:,:,1) = tmask(:,:,1)     ! surface 
     287      wumask(:,:,1) = umask(:,:,1) 
     288      wvmask(:,:,1) = vmask(:,:,1) 
     289      DO jk = 2, jpk                   ! interior values 
     290         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 
     291         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)    
     292         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 
    293293      END DO 
    294294 
     
    377377      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration 
    378378         !                                                 ! Increased lateral friction near of some straits 
    379          IF( nn_cla == 0 ) THEN 
    380             !                                ! Gibraltar strait  : partial slip (fmask=0.5) 
    381             ij0 = 101   ;   ij1 = 101 
    382             ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
    383             ij0 = 102   ;   ij1 = 102 
    384             ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
    385             ! 
    386             !                                ! Bab el Mandeb : partial slip (fmask=1) 
    387             ij0 =  87   ;   ij1 =  88 
    388             ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
    389             ij0 =  88   ;   ij1 =  88 
    390             ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
    391             ! 
    392          ENDIF 
     379         !                                ! Gibraltar strait  : partial slip (fmask=0.5) 
     380         ij0 = 101   ;   ij1 = 101 
     381         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
     382         ij0 = 102   ;   ij1 = 102 
     383         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
     384         ! 
     385         !                                ! Bab el Mandeb : partial slip (fmask=1) 
     386         ij0 =  87   ;   ij1 =  88 
     387         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
     388         ij0 =  88   ;   ij1 =  88 
     389         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
     390         ! 
    393391         !                                ! Danish straits  : strong slip (fmask > 2) 
    394392! We keep this as an example but it is instable in this case  
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90

    r3294 r5770  
    3737      !!                -> not good if located at too high latitude... 
    3838      !!---------------------------------------------------------------------- 
    39       ! 
    4039      REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point 
    4140      INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point 
     
    4948      IF( nn_timing == 1 )  CALL timing_start('dom_ngb') 
    5049      ! 
    51       CALL wrk_alloc( jpi, jpj, zglam, zgphi, zmask, zdist ) 
     50      CALL wrk_alloc( jpi,jpj,  zglam, zgphi, zmask, zdist ) 
    5251      ! 
    5352      zmask(:,:) = 0._wp 
     
    7675      ENDIF 
    7776      ! 
    78       CALL wrk_dealloc( jpi, jpj, zglam, zgphi, zmask, zdist ) 
     77      CALL wrk_dealloc( jpi,jpj,  zglam, zgphi, zmask, zdist ) 
    7978      ! 
    8079      IF( nn_timing == 1 )  CALL timing_stop('dom_ngb') 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r5737 r5770  
    3636   !!---------------------------------------------------------------------- 
    3737CONTAINS 
    38  
    39  
    4038 
    4139   SUBROUTINE dom_wri_coordinate 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5656 r5770  
    501501            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    502502               !                                             ! ===================== 
    503                IF( nn_cla == 0 ) THEN 
    504                   ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open  
    505                   ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995) 
    506                   DO ji = mi0(ii0), mi1(ii1) 
    507                      DO jj = mj0(ij0), mj1(ij1) 
    508                         mbathy(ji,jj) = 15 
    509                      END DO 
     503               ! 
     504               ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open  
     505               ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995) 
     506               DO ji = mi0(ii0), mi1(ii1) 
     507                  DO jj = mj0(ij0), mj1(ij1) 
     508                     mbathy(ji,jj) = 15 
    510509                  END DO 
    511                   IF(lwp) WRITE(numout,*) 
    512                   IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
    513                   ! 
    514                   ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open 
    515                   ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995) 
    516                   DO ji = mi0(ii0), mi1(ii1) 
    517                      DO jj = mj0(ij0), mj1(ij1) 
    518                         mbathy(ji,jj) = 12 
    519                      END DO 
     510               END DO 
     511               IF(lwp) WRITE(numout,*) 
     512               IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
     513               ! 
     514               ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open 
     515               ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995) 
     516               DO ji = mi0(ii0), mi1(ii1) 
     517                  DO jj = mj0(ij0), mj1(ij1) 
     518                     mbathy(ji,jj) = 12 
    520519                  END DO 
    521                   IF(lwp) WRITE(numout,*) 
    522                   IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    523                ENDIF 
     520               END DO 
     521               IF(lwp) WRITE(numout,*) 
     522               IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    524523               ! 
    525524            ENDIF 
     
    545544            !        
    546545            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    547                ! 
    548               IF( nn_cla == 0 ) THEN 
    549                  ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open  
    550                  ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995) 
    551                  DO ji = mi0(ii0), mi1(ii1) 
    552                     DO jj = mj0(ij0), mj1(ij1) 
    553                        bathy(ji,jj) = 284._wp 
    554                     END DO 
     546            ! 
     547              ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open  
     548              ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995) 
     549              DO ji = mi0(ii0), mi1(ii1) 
     550                 DO jj = mj0(ij0), mj1(ij1) 
     551                    bathy(ji,jj) = 284._wp 
    555552                 END DO 
    556                  IF(lwp) WRITE(numout,*)      
    557                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
    558                  ! 
    559                  ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open 
    560                  ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995) 
    561                  DO ji = mi0(ii0), mi1(ii1) 
    562                     DO jj = mj0(ij0), mj1(ij1) 
    563                        bathy(ji,jj) = 137._wp 
    564                     END DO 
     553               END DO 
     554              IF(lwp) WRITE(numout,*)      
     555              IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 
     556              ! 
     557              ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open 
     558              ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995) 
     559               DO ji = mi0(ii0), mi1(ii1) 
     560                 DO jj = mj0(ij0), mj1(ij1) 
     561                    bathy(ji,jj) = 137._wp 
    565562                 END DO 
    566                  IF(lwp) WRITE(numout,*) 
    567                  IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    568               ENDIF 
     563              END DO 
     564              IF(lwp) WRITE(numout,*) 
     565               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 
    569566              ! 
    570567           ENDIF 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r5215 r5770  
    174174            END DO 
    175175         END DO 
    176          IF( nn_cla == 1 ) THEN                          ! Cross Land advection 
    177             il0 = 138   ;   il1 = 138                          ! set T & S profile at Gibraltar Strait 
    178             ij0 = 101   ;   ij1 = 102 
    179             ii0 = 139   ;   ii1 = 139 
    180             DO jl = mi0(il0), mi1(il1) 
    181                DO jj = mj0(ij0), mj1(ij1) 
    182                   DO ji = mi0(ii0), mi1(ii1) 
    183                      sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:) 
    184                      sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:) 
    185                   END DO 
    186                END DO 
    187             END DO 
    188             il0 = 164   ;   il1 = 164                          ! set T & S profile at Bab el Mandeb Strait 
    189             ij0 =  87   ;   ij1 =  88 
    190             ii0 = 161   ;   ii1 = 163 
    191             DO jl = mi0(il0), mi1(il1) 
    192                DO jj = mj0(ij0), mj1(ij1) 
    193                   DO ji = mi0(ii0), mi1(ii1) 
    194                      sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:) 
    195                      sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:) 
    196                   END DO 
    197                END DO 
    198             END DO 
    199          ELSE                                            ! No Cross Land advection 
    200             ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea 
    201             ii0 = 148   ;   ii1 = 160 
    202             sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp 
    203             sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp 
    204             sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp 
    205          ENDIF 
     176         ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea 
     177         ii0 = 148   ;   ii1 = 160 
     178         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp 
     179         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp 
     180         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp 
    206181      ENDIF 
    207182      ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5766 r5770  
    3434   USE dtatsd          ! data temperature and salinity   (dta_tsd routine) 
    3535   USE dtauvd          ! data: U & V current             (dta_uvd routine) 
    36    USE zpshde          ! partial step: hor. derivative (zps_hde routine) 
    37    USE eosbn2          ! equation of state            (eos bn2 routine) 
    3836   USE domvvl          ! varying vertical mesh 
    3937   USE dynspg_oce      ! pressure gradient schemes 
     
    119117            ENDIF 
    120118            IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 
    121                CALL wrk_alloc( jpi, jpj, jpk, 2, zuvd ) 
     119               CALL wrk_alloc( jpi,jpj,jpk,2,  zuvd ) 
    122120               CALL dta_uvd( nit000, zuvd ) 
    123121               ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:) 
    124122               vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:) 
    125                CALL wrk_dealloc( jpi, jpj, jpk, 2, zuvd ) 
     123               CALL wrk_dealloc( jpi,jpj,jpk,2,  zuvd ) 
    126124            ENDIF 
    127125         ENDIF 
    128          ! 
     126         !    
    129127         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 
    130128         IF( lk_vvl ) THEN 
    131129            DO jk = 1, jpk 
    132130               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
    133             ENDDO 
     131            END DO 
    134132         ENDIF 
    135133         !  
     
    456454      !!---------------------------------------------------------------------- 
    457455      ! 
    458       CALL wrk_alloc( jpi, jpj, jpk, zprn) 
     456      CALL wrk_alloc( jpi,jpj,jpk,  zprn) 
    459457      ! 
    460458      IF(lwp) WRITE(numout,*)  
     
    553551      rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value 
    554552      ! 
    555       CALL wrk_dealloc( jpi, jpj, jpk, zprn) 
     553      CALL wrk_dealloc( jpi,jpj,jpk,  zprn) 
    556554      ! 
    557555   END SUBROUTINE istate_uvg 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r5737 r5770  
    1818   !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here 
    1919   !!            3.6  ! 2014-11  (P. Mathiot)          isf            added directly here 
     20   !!            3.7  ! 2015-10  (G. Madec) remove cross-land advection 
    2021   !!---------------------------------------------------------------------- 
    2122 
     
    2930   USE sbcrnf          ! river runoff  
    3031   USE sbcisf          ! ice shelf  
    31    USE cla             ! cross land advection             (cla_div routine) 
    3232   USE in_out_manager  ! I/O manager 
    3333   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    6868      !!         - compute the now divergence given by : 
    6969      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    70       !!      correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf) 
    71       !!      and cross land flow (div_cla)  
     70      !!      correct hdiv with runoff inflow (div_rnf) and ice shelf melting (div_isf) 
    7271      !!              II. vorticity : 
    7372      !!         - save the curl computed at the previous time-step 
     
    226225      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field) 
    227226      IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
    228       IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence) 
    229227       
    230228      ! 4. Lateral boundary conditions on hdivn and rotn 
     
    256254      !!      - compute the now divergence given by : 
    257255      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    258       !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)  
     256      !!      correct hdiv with runoff inflow (div_rnf) 
    259257      !!              - Relavtive Vorticity : 
    260258      !!      - save the curl computed at the previous time-step (rotb = rotn) 
     
    325323      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field) 
    326324      IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
    327       IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    328325      ! 
    329326      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change) 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r5120 r5770  
    271271      CALL solver_init( nit000 )   ! Elliptic solver initialisation 
    272272#endif 
    273  
    274       !                        ! Control of timestep choice 
    275       IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN 
    276          IF( nn_cla == 1 )   CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' ) 
    277       ENDIF 
    278  
    279273      !               ! Control of hydrostatic pressure choice 
    280       IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN 
    281          CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 
    282       ENDIF 
     274      IF( lk_dynspg_ts .AND. ln_dynhpg_imp )   CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 
    283275      ! 
    284276      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init') 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r4990 r5770  
    1414   !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 
    1515   !!            3.7  !  2014-04  (F. Roquet, G. Madec)  add some trends diag 
     16   !!             -   !  2014-12  (G. Madec) remove cross-land advection (cla) 
    1617   !!---------------------------------------------------------------------- 
    1718#if defined key_dynspg_flt   ||   defined key_esopa   
     
    3637   USE bdydyn          ! ocean open boundary condition on dynamics 
    3738   USE bdyvol          ! ocean open boundary condition (bdy_vol routine) 
    38    USE cla             ! cross land advection 
    3939   USE trd_oce         ! trends: ocean variables 
    4040   USE trddyn          ! trend manager: dynamics 
     
    206206      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces  
    207207#endif 
    208       IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va)) 
    209208 
    210209      ! compute the next vertically averaged velocity (effect of the additional force not included) 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r5737 r5770  
    8080      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
    8181      ! 
    82       CALL wrk_alloc( jpi, jpj, zhdiv )  
     82      CALL wrk_alloc( jpi,jpj,  zhdiv )  
    8383      ! 
    8484      IF( kt == nit000 ) THEN 
    85          ! 
    8685         IF(lwp) WRITE(numout,*) 
    8786         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 
    8887         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    89          ! 
    9088      ENDIF 
    9189      ! 
     
    159157      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    160158      !!---------------------------------------------------------------------- 
    161       ! 
    162       INTEGER, INTENT(in) ::   kt           ! time step 
     159      INTEGER, INTENT(in) ::   kt   ! time step 
     160      ! 
     161      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     162      REAL(wp) ::   z1_2dt       ! local scalars 
    163163      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
    164164      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv 
    165       ! 
    166       INTEGER             ::   ji, jj, jk   ! dummy loop indices 
    167       REAL(wp)            ::   z1_2dt       ! local scalars 
    168       !!---------------------------------------------------------------------- 
    169        
     165      !!---------------------------------------------------------------------- 
     166      ! 
    170167      IF( nn_timing == 1 )  CALL timing_start('wzv') 
    171168      ! 
    172169      IF( kt == nit000 ) THEN 
    173          ! 
    174170         IF(lwp) WRITE(numout,*) 
    175171         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 
     
    177173         ! 
    178174         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
    179          ! 
    180175      ENDIF 
    181176      !                                           !------------------------------! 
     
    216211 
    217212#if defined key_bdy 
    218       IF (lk_bdy) THEN 
     213      IF( lk_bdy ) THEN 
    219214         DO jk = 1, jpkm1 
    220215            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
     
    224219      ! 
    225220      IF( nn_timing == 1 )  CALL timing_stop('wzv') 
    226  
    227  
     221      ! 
    228222   END SUBROUTINE wzv 
     223 
    229224 
    230225   SUBROUTINE ssh_swp( kt ) 
     
    265260         sshb(:,:) = sshn(:,:)                           ! before <-- now 
    266261         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now) 
     262         ! 
    267263      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap 
    268264         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r5758 r5770  
    77   !!            3.3  !  2010-09  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport 
    88   !!            3.6  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation 
    9    !!            4.0  !  2014-05  (G. Madec)  Add 2nd/4th order cases for CEN and FCT schemes  
     9   !!            3.7  !  2014-05  (G. Madec)  Add 2nd/4th order cases for CEN and FCT schemes  
    1010   !!             -   !  2014-12  (G. Madec) suppression of cross land advection option 
    1111   !!---------------------------------------------------------------------- 
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   tra_adv      : compute ocean tracer advection trend 
    15    !!   tra_adv_ctl  : control the different options of advection scheme 
    16    !!---------------------------------------------------------------------- 
    17    USE oce             ! ocean dynamics and active tracers 
    18    USE dom_oce         ! ocean space and time domain 
    19    USE domvvl          ! variable vertical scale factors 
    20    USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine) 
    21    USE traadv_tvd      ! TVD      scheme           (tra_adv_tvd    routine) 
    22    USE traadv_muscl    ! MUSCL    scheme           (tra_adv_muscl  routine) 
    23    USE traadv_muscl2   ! MUSCL2   scheme           (tra_adv_muscl2 routine) 
    24    USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine) 
    25    USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine) 
    26    USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine) 
    27    USE cla             ! cross land advection      (cla_traadv     routine) 
    28    USE ldftra          ! lateral diffusion coefficient on tracers 
    29    USE ldfslp          ! Lateral diffusion: slopes of neutral surfaces 
     14   !!   tra_adv       : compute ocean tracer advection trend 
     15   !!   tra_adv_ctl   : control the different options of advection scheme 
     16   !!---------------------------------------------------------------------- 
     17   USE oce            ! ocean dynamics and active tracers 
     18   USE dom_oce        ! ocean space and time domain 
     19   USE domvvl         ! variable vertical scale factors 
     20   USE traadv_cen     ! centered scheme           (tra_adv_cen  routine) 
     21   USE traadv_fct     ! FCT      scheme           (tra_adv_fct  routine) 
     22   USE traadv_mus     ! MUSCL    scheme           (tra_adv_mus  routine) 
     23   USE traadv_ubs     ! UBS      scheme           (tra_adv_ubs  routine) 
     24   USE traadv_qck     ! QUICKEST scheme           (tra_adv_qck  routine) 
     25   USE traadv_mle     ! ML eddy induced velocity  (tra_adv_mle  routine) 
     26   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff. 
     27   USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces 
    3028   ! 
    31    USE in_out_manager  ! I/O manager 
    32    USE iom             ! I/O module 
    33    USE prtctl          ! Print control 
    34    USE lib_mpp         ! MPP library 
    35    USE wrk_nemo        ! Memory Allocation 
    36    USE timing          ! Timing 
    37    USE sbc_oce 
     29   USE in_out_manager ! I/O manager 
     30   USE iom            ! I/O module 
     31   USE prtctl         ! Print control 
     32   USE lib_mpp        ! MPP library 
     33   USE wrk_nemo       ! Memory Allocation 
     34   USE timing         ! Timing 
     35 
    3836   USE diaptr          ! Poleward heat transport  
    39  
    4037 
    4138   IMPLICIT NONE 
     
    4542   PUBLIC   tra_adv_init   ! routine called by opa module 
    4643 
    47    !                              !!* Namelist namtra_adv * 
    48    LOGICAL ::   ln_traadv_cen2     ! 2nd order centered scheme flag 
    49    LOGICAL ::   ln_traadv_tvd      ! TVD scheme flag 
    50    LOGICAL ::   ln_traadv_tvd_zts  ! TVD scheme flag with vertical sub time-stepping 
    51    LOGICAL ::   ln_traadv_muscl    ! MUSCL scheme flag 
    52    LOGICAL ::   ln_traadv_muscl2   ! MUSCL2 scheme flag 
    53    LOGICAL ::   ln_traadv_ubs      ! UBS scheme flag 
    54    LOGICAL ::   ln_traadv_qck      ! QUICKEST scheme flag 
    55    LOGICAL ::   ln_traadv_msc_ups  ! use upstream scheme within muscl 
    56  
    57  
    58    INTEGER ::   nadv   ! choice of the type of advection scheme 
    59  
     44   !                            !!* Namelist namtra_adv * 
     45   LOGICAL ::   ln_traadv_cen    ! centered scheme flag 
     46   INTEGER ::      nn_cen_h, nn_cen_v   ! =2/4 : horizontal and vertical choices of the order of CEN scheme 
     47   LOGICAL ::   ln_traadv_fct    ! FCT scheme flag 
     48   INTEGER ::      nn_fct_h, nn_fct_v   ! =2/4 : horizontal and vertical choices of the order of FCT scheme 
     49   INTEGER ::      nn_fct_zts           ! >=1 : 2nd order FCT with vertical sub-timestepping 
     50   LOGICAL ::   ln_traadv_mus    ! MUSCL scheme flag 
     51   LOGICAL ::      ln_mus_ups           ! use upstream scheme in vivcinity of river mouths 
     52   LOGICAL ::   ln_traadv_ubs    ! UBS scheme flag 
     53   INTEGER ::      nn_ubs_v             ! =2/4 : vertical choice of the order of UBS scheme 
     54   LOGICAL ::   ln_traadv_qck    ! QUICKEST scheme flag 
     55 
     56   INTEGER ::              nadv             ! choice of the type of advection scheme 
     57   ! 
     58   !                                        ! associated indices: 
     59   INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection 
     60   INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme 
     61   INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme 
     62   INTEGER, PARAMETER ::   np_FCT_zts = 3   ! 2nd order FCT scheme with vertical sub-timestepping 
     63   INTEGER, PARAMETER ::   np_MUS     = 4   ! MUSCL scheme 
     64   INTEGER, PARAMETER ::   np_UBS     = 5   ! 3rd order Upstream Biased Scheme 
     65   INTEGER, PARAMETER ::   np_QCK     = 6   ! QUICK scheme 
     66    
    6067   !! * Substitutions 
    6168#  include "domzgr_substitute.h90" 
    6269#  include "vectopt_loop_substitute.h90" 
    6370   !!---------------------------------------------------------------------- 
    64    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     71   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    6572   !! $Id$ 
    6673   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    8491      IF( nn_timing == 1 )  CALL timing_start('tra_adv') 
    8592      ! 
    86       CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn ) 
     93      CALL wrk_alloc( jpi,jpj,jpk,  zun, zvn, zwn ) 
    8794      ! 
    8895      !                                          ! set time step 
     
    93100      ENDIF 
    94101      ! 
    95       IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_traadv( kt )       !==  Cross Land Advection  ==! (hor. advection) 
    96       ! 
    97       !                                               !==  effective transport  ==! 
     102      !                                         !==  effective transport  ==! 
    98103      DO jk = 1, jpkm1 
    99104         zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only 
     
    102107      END DO 
    103108      ! 
    104       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     109      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections 
    105110         zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 
    106111         zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 
    107112      ENDIF 
    108113      ! 
    109       zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    110       zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    111       zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
     114      zun(:,:,jpk) = 0._wp                                                      ! no transport trough the bottom 
     115      zvn(:,:,jpk) = 0._wp 
     116      zwn(:,:,jpk) = 0._wp 
    112117      ! 
    113118      IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   & 
     
    120125      CALL iom_put( "wocetr_eff", zwn ) 
    121126      ! 
     127!!gm ??? 
    122128      IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF  
    123       ! 
    124     
    125       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==! 
    126       CASE ( 1 )   ;    CALL tra_adv_cen2   ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered 
    127       CASE ( 2 )   ;    CALL tra_adv_tvd    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD  
    128       CASE ( 3 )   ;    CALL tra_adv_muscl  ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts, ln_traadv_msc_ups )   !  MUSCL  
    129       CASE ( 4 )   ;    CALL tra_adv_muscl2 ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2  
    130       CASE ( 5 )   ;    CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS  
    131       CASE ( 6 )   ;    CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST  
    132       CASE ( 7 )   ;    CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD ZTS 
    133       ! 
    134       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==! 
    135          CALL tra_adv_cen2  ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    136          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv0 - Ta: ', mask1=tmask,               & 
    137             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    138          CALL tra_adv_tvd   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    139          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv1 - Ta: ', mask1=tmask,               & 
    140             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    141          CALL tra_adv_muscl ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts, ln_traadv_msc_ups )           
    142          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv3 - Ta: ', mask1=tmask,               & 
    143             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    144          CALL tra_adv_muscl2( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    145          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv4 - Ta: ', mask1=tmask,               & 
    146             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    147          CALL tra_adv_ubs   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    148          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv5 - Ta: ', mask1=tmask,               & 
    149             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    150          CALL tra_adv_qck   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    151          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv6 - Ta: ', mask1=tmask,               & 
    152             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     129!!gm ??? 
     130      ! 
     131      SELECT CASE ( nadv )                      !==  compute advection trend and add it to general trend  ==! 
     132      ! 
     133      CASE ( np_CEN )                                    ! Centered scheme : 2nd / 4th order 
     134         CALL tra_adv_cen    ( kt, nit000, 'TRA',         zun, zvn, zwn     , tsn, tsa, jpts, nn_cen_h, nn_cen_v ) 
     135      CASE ( np_FCT )                                    ! FCT scheme      : 2nd / 4th order 
     136         CALL tra_adv_fct    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v ) 
     137      CASE ( np_FCT_zts )                                ! 2nd order FCT with vertical time-splitting 
     138         CALL tra_adv_fct_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_fct_zts ) 
     139      CASE ( np_MUS )                                    ! MUSCL 
     140         CALL tra_adv_mus    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts        , ln_mus_ups )  
     141      CASE ( np_UBS )                                    ! UBS 
     142         CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_ubs_v   ) 
     143      CASE ( np_QCK )                                    ! QUICKEST 
     144         CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts                     ) 
     145      ! 
    153146      END SELECT 
    154147      ! 
     
    159152      IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv' ) 
    160153      ! 
    161       CALL wrk_dealloc( jpi, jpj, jpk, zun, zvn, zwn ) 
     154      CALL wrk_dealloc( jpi,jpj,jpk,  zun, zvn, zwn ) 
    162155      !                                           
    163156   END SUBROUTINE tra_adv 
     
    171164      !!              tracer advection schemes and set nadv 
    172165      !!---------------------------------------------------------------------- 
    173       INTEGER ::   ioptio 
    174       INTEGER ::   ios                 ! Local integer output status for namelist read 
    175       !! 
    176       NAMELIST/namtra_adv/ ln_traadv_cen2 , ln_traadv_tvd,     & 
    177          &                 ln_traadv_muscl, ln_traadv_muscl2,  & 
    178          &                 ln_traadv_ubs  , ln_traadv_qck,     & 
    179          &                 ln_traadv_msc_ups, ln_traadv_tvd_zts 
    180       !!---------------------------------------------------------------------- 
    181  
    182       REWIND( numnam_ref )              ! Namelist namtra_adv in reference namelist : Tracer advection scheme 
     166      INTEGER ::   ioptio, ios   ! Local integers 
     167      ! 
     168      NAMELIST/namtra_adv/ ln_traadv_cen, nn_cen_h, nn_cen_v,               &   ! CEN 
     169         &                 ln_traadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts,   &   ! FCT 
     170         &                 ln_traadv_mus,                     ln_mus_ups,   &   ! MUSCL 
     171         &                 ln_traadv_ubs,           nn_ubs_v,               &   ! UBS 
     172         &                 ln_traadv_qck                                        ! QCK 
     173      !!---------------------------------------------------------------------- 
     174      ! 
     175      !                                !==  Namelist  ==! 
     176      ! 
     177      REWIND( numnam_ref )                   ! Namelist namtra_adv in reference namelist : Tracer advection scheme 
    183178      READ  ( numnam_ref, namtra_adv, IOSTAT = ios, ERR = 901) 
    184179901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in reference namelist', lwp ) 
    185  
    186       REWIND( numnam_cfg )              ! Namelist namtra_adv in configuration namelist : Tracer advection scheme 
     180      ! 
     181      REWIND( numnam_cfg )                   ! Namelist namtra_adv in configuration namelist : Tracer advection scheme 
    187182      READ  ( numnam_cfg, namtra_adv, IOSTAT = ios, ERR = 902 ) 
    188183902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in configuration namelist', lwp ) 
    189184      IF(lwm) WRITE ( numond, namtra_adv ) 
    190185 
    191       IF(lwp) THEN                    ! Namelist print 
     186      IF(lwp) THEN                           ! Namelist print 
    192187         WRITE(numout,*) 
    193188         WRITE(numout,*) 'tra_adv_init : choice/control of the tracer advection scheme' 
    194189         WRITE(numout,*) '~~~~~~~~~~~' 
    195190         WRITE(numout,*) '   Namelist namtra_adv : chose a advection scheme for tracers' 
    196          WRITE(numout,*) '      2nd order advection scheme     ln_traadv_cen2    = ', ln_traadv_cen2 
    197          WRITE(numout,*) '      TVD advection scheme           ln_traadv_tvd     = ', ln_traadv_tvd 
    198          WRITE(numout,*) '      MUSCL  advection scheme        ln_traadv_muscl   = ', ln_traadv_muscl 
    199          WRITE(numout,*) '      MUSCL2 advection scheme        ln_traadv_muscl2  = ', ln_traadv_muscl2 
    200          WRITE(numout,*) '      UBS    advection scheme        ln_traadv_ubs     = ', ln_traadv_ubs 
    201          WRITE(numout,*) '      QUICKEST advection scheme      ln_traadv_qck     = ', ln_traadv_qck 
    202          WRITE(numout,*) '      upstream scheme within muscl   ln_traadv_msc_ups = ', ln_traadv_msc_ups 
    203          WRITE(numout,*) '      TVD advection scheme with zts  ln_traadv_tvd_zts = ', ln_traadv_tvd_zts 
    204       ENDIF 
    205  
    206       ioptio = 0                      ! Parameter control 
    207       IF( ln_traadv_cen2   )   ioptio = ioptio + 1 
    208       IF( ln_traadv_tvd    )   ioptio = ioptio + 1 
    209       IF( ln_traadv_muscl  )   ioptio = ioptio + 1 
    210       IF( ln_traadv_muscl2 )   ioptio = ioptio + 1 
    211       IF( ln_traadv_ubs    )   ioptio = ioptio + 1 
    212       IF( ln_traadv_qck    )   ioptio = ioptio + 1 
    213       IF( ln_traadv_tvd_zts)   ioptio = ioptio + 1 
    214       IF( lk_esopa         )   ioptio =          1 
    215  
    216       IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts )   & 
    217          .AND. ln_isfcav )   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity') 
    218  
    219       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtra_adv' ) 
    220  
    221       !                              ! Set nadv 
    222       IF( ln_traadv_cen2   )   nadv =  1 
    223       IF( ln_traadv_tvd    )   nadv =  2 
    224       IF( ln_traadv_muscl  )   nadv =  3 
    225       IF( ln_traadv_muscl2 )   nadv =  4 
    226       IF( ln_traadv_ubs    )   nadv =  5 
    227       IF( ln_traadv_qck    )   nadv =  6 
    228       IF( ln_traadv_tvd_zts)   nadv =  7 
    229       IF( lk_esopa         )   nadv = -1 
    230  
    231       IF(lwp) THEN                   ! Print the choice 
     191         WRITE(numout,*) '      centered scheme                           ln_traadv_cen = ', ln_traadv_cen 
     192         WRITE(numout,*) '            horizontal 2nd/4th order               nn_cen_h   = ', nn_fct_h 
     193         WRITE(numout,*) '            vertical   2nd/4th order               nn_cen_v   = ', nn_fct_v 
     194         WRITE(numout,*) '      Flux Corrected Transport scheme           ln_traadv_fct = ', ln_traadv_fct 
     195         WRITE(numout,*) '            horizontal 2nd/4th order               nn_fct_h   = ', nn_fct_h 
     196         WRITE(numout,*) '            vertical   2nd/4th order               nn_fct_v   = ', nn_fct_v 
     197         WRITE(numout,*) '            2nd order + vertical sub-timestepping  nn_fct_zts = ', nn_fct_zts 
     198         WRITE(numout,*) '      MUSCL scheme                              ln_traadv_mus = ', ln_traadv_mus 
     199         WRITE(numout,*) '            + upstream scheme near river mouths    ln_mus_ups = ', ln_mus_ups 
     200         WRITE(numout,*) '      UBS scheme                                ln_traadv_ubs = ', ln_traadv_ubs 
     201         WRITE(numout,*) '            vertical   2nd/4th order               nn_ubs_v   = ', nn_ubs_v 
     202         WRITE(numout,*) '      QUICKEST scheme                           ln_traadv_qck = ', ln_traadv_qck 
     203      ENDIF 
     204 
     205      ioptio = 0                       !==  Parameter control  ==! 
     206      IF( ln_traadv_cen )   ioptio = ioptio + 1 
     207      IF( ln_traadv_fct )   ioptio = ioptio + 1 
     208      IF( ln_traadv_mus )   ioptio = ioptio + 1 
     209      IF( ln_traadv_ubs )   ioptio = ioptio + 1 
     210      IF( ln_traadv_qck )   ioptio = ioptio + 1 
     211      ! 
     212      IF( ioptio == 0 ) THEN 
     213         nadv = np_NO_adv 
     214         CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' ) 
     215      ENDIF 
     216      IF( ioptio /= 1 )   CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 
     217      ! 
     218      IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 )   &          ! Centered 
     219                        .AND. ( nn_cen_v /= 2 .AND. nn_cen_v /= 4 )   ) THEN 
     220        CALL ctl_stop( 'tra_adv_init: CEN scheme, choose 2nd or 4th order' ) 
     221      ENDIF 
     222      IF( ln_traadv_fct .AND. ( nn_fct_h /= 2 .AND. nn_fct_h /= 4 )   &          ! FCT 
     223                        .AND. ( nn_fct_v /= 2 .AND. nn_fct_v /= 4 )   ) THEN 
     224        CALL ctl_stop( 'tra_adv_init: FCT scheme, choose 2nd or 4th order' ) 
     225      ENDIF 
     226      IF( ln_traadv_fct .AND. nn_fct_zts > 0 ) THEN 
     227         IF( nn_fct_h == 4 ) THEN 
     228            nn_fct_h = 2 
     229            CALL ctl_stop( 'tra_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' ) 
     230         ENDIF 
     231         IF( lk_vvl ) THEN 
     232            CALL ctl_stop( 'tra_adv_init: vertical sub-timestepping not allow in non-linear free surface' ) 
     233         ENDIF 
     234         IF( nn_fct_zts == 1 )   CALL ctl_warn( 'tra_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' ) 
     235      ENDIF 
     236      IF( ln_traadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 )   ) THEN     ! UBS 
     237        CALL ctl_stop( 'tra_adv_init: UBS scheme, choose 2nd or 4th order' ) 
     238      ENDIF 
     239      IF( ln_traadv_ubs .AND. nn_ubs_v == 4 ) THEN 
     240         CALL ctl_warn( 'tra_adv_init: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' ) 
     241      ENDIF 
     242      IF( ln_isfcav ) THEN                                                       ! ice-shelf cavities 
     243         IF(  ln_traadv_cen .AND. nn_cen_v /= 4    .OR.   &                            ! NO 4th order with ISF 
     244            & ln_traadv_fct .AND. nn_fct_v /= 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 
     245      ENDIF 
     246      ! 
     247      !                                !==  used advection scheme  ==!   
     248      !                                      ! set nadv 
     249      IF( ln_traadv_cen                      )   nadv = np_CEN 
     250      IF( ln_traadv_fct                      )   nadv = np_FCT 
     251      IF( ln_traadv_fct .AND. nn_fct_zts > 0 )   nadv = np_FCT_zts 
     252      IF( ln_traadv_mus                      )   nadv = np_MUS 
     253      IF( ln_traadv_ubs                      )   nadv = np_UBS 
     254      IF( ln_traadv_qck                      )   nadv = np_QCK 
     255 
     256      IF(lwp) THEN                           ! Print the choice 
    232257         WRITE(numout,*) 
    233          IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used' 
    234          IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used' 
    235          IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used' 
    236          IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    scheme is used' 
    237          IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used' 
    238          IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used' 
    239          IF( nadv ==  7 )   WRITE(numout,*) '         TVD ZTS   scheme is used' 
    240          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme' 
    241       ENDIF 
    242       ! 
    243       CALL tra_adv_mle_init          ! initialisation of the Mixed Layer Eddy parametrisation (MLE) 
     258         IF( nadv == np_NO_adv  )   WRITE(numout,*) '         NO T-S advection' 
     259         IF( nadv == np_CEN     )   WRITE(numout,*) '         CEN      scheme is used. Horizontal order: ', nn_cen_h,   & 
     260            &                                                                        ' Vertical   order: ', nn_cen_v 
     261         IF( nadv == np_FCT     )   WRITE(numout,*) '         FCT      scheme is used. Horizontal order: ', nn_fct_h,   & 
     262            &                                                                        ' Vertical   order: ', nn_fct_v 
     263         IF( nadv == np_FCT_zts )   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 
     264         IF( nadv == np_MUS     )   WRITE(numout,*) '         MUSCL    scheme is used' 
     265         IF( nadv == np_UBS     )   WRITE(numout,*) '         UBS      scheme is used' 
     266         IF( nadv == np_QCK     )   WRITE(numout,*) '         QUICKEST scheme is used' 
     267      ENDIF 
     268      ! 
     269      CALL tra_adv_mle_init            !== initialisation of the Mixed Layer Eddy parametrisation (MLE)  ==! 
    244270      ! 
    245271   END SUBROUTINE tra_adv_init 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

    r5737 r5770  
    1 MODULE traadv_tvd 
     1MODULE traadv_fct 
    22   !!============================================================================== 
    3    !!                       ***  MODULE  traadv_tvd  *** 
    4    !! Ocean  tracers:  horizontal & vertical advective trend 
     3   !!                       ***  MODULE  traadv_fct  *** 
     4   !! Ocean  tracers:  horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method) 
    55   !!============================================================================== 
    6    !! History :  OPA  !  1995-12  (L. Mortier)  Original code 
    7    !!                 !  2000-01  (H. Loukos)  adapted to ORCA  
    8    !!                 !  2000-10  (MA Foujols E.Kestenare)  include file not routine 
    9    !!                 !  2000-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes 
    10    !!                 !  2001-07  (E. Durand G. Madec)  adaptation to ORCA config 
    11    !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module 
    12    !!    NEMO    1.0  !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl 
    13    !!            2.0  !  2008-04  (S. Cravatte) add the i-, j- & k- trends computation 
    14    !!             -   !  2009-11  (V. Garnier) Surface pressure gradient organization 
    15    !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport 
     6   !! History :  3.7  !  2015-09  (L. Debreu, G. Madec)  original code (inspired from traadv_tvd.F90) 
    167   !!---------------------------------------------------------------------- 
    178 
    189   !!---------------------------------------------------------------------- 
    19    !!   tra_adv_tvd   : update the tracer trend with the 3D advection trends using a TVD scheme 
    20    !!   nonosc        : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     10   !!  tra_adv_fct    : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme 
     11   !!  tra_adv_fct_zts: update the tracer trend with a 3D advective trends using a 2nd order FCT scheme  
     12   !!                   with sub-time-stepping in the vertical direction 
     13   !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     14   !!  interp_4th_cpt : 4th order compact scheme for the vertical component of the advection 
    2115   !!---------------------------------------------------------------------- 
    2216   USE oce            ! ocean dynamics and active tracers 
     
    2822   USE diaptr         ! poleward transport diagnostics 
    2923   ! 
     24   USE in_out_manager ! I/O manager 
    3025   USE lib_mpp        ! MPP library 
    3126   USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    32    USE in_out_manager ! I/O manager 
     27   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3328   USE wrk_nemo       ! Memory Allocation 
    3429   USE timing         ! Timing 
    35    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3630 
    3731   IMPLICIT NONE 
    3832   PRIVATE 
    3933 
    40    PUBLIC   tra_adv_tvd        ! routine called by traadv.F90 
    41    PUBLIC   tra_adv_tvd_zts    ! routine called by traadv.F90 
    42  
    43    LOGICAL ::   l_trd   ! flag to compute trends 
     34   PUBLIC   tra_adv_fct        ! routine called by traadv.F90 
     35   PUBLIC   tra_adv_fct_zts    ! routine called by traadv.F90 
     36   PUBLIC   interp_4th_cpt     ! routine called by traadv_cen.F90 
     37 
     38   LOGICAL  ::   l_trd   ! flag to compute trends 
     39   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6 
    4440 
    4541   !! * Substitutions 
     
    4743#  include "vectopt_loop_substitute.h90" 
    4844   !!---------------------------------------------------------------------- 
    49    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     45   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    5046   !! $Id$ 
    5147   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    5349CONTAINS 
    5450 
    55    SUBROUTINE tra_adv_tvd ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    56       &                                       ptb, ptn, pta, kjpt ) 
    57       !!---------------------------------------------------------------------- 
    58       !!                  ***  ROUTINE tra_adv_tvd  *** 
     51   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn,       & 
     52      &                                              ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v ) 
     53      !!---------------------------------------------------------------------- 
     54      !!                  ***  ROUTINE tra_adv_fct  *** 
    5955      !!  
    6056      !! **  Purpose :   Compute the now trend due to total advection of  
    6157      !!       tracers and add it to the general trend of tracer equations 
    6258      !! 
    63       !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with 
    64       !!       corrected flux (monotonic correction) 
    65       !!       note: - this advection scheme needs a leap-frog time scheme 
     59      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction 
     60      !!               (choice through the value of kn_fct) 
     61      !!               - 4th order compact scheme on the vertical  
     62      !!               - corrected flux (monotonic correction)  
    6663      !! 
    6764      !! ** Action : - update (pta) with the now advective tracer trends 
    68       !!             - save the trends  
    69       !!---------------------------------------------------------------------- 
    70       USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace 
    71       ! 
     65      !!             - send the trends for further diagnostics 
     66      !!---------------------------------------------------------------------- 
    7267      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    7368      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    7469      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    7570      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
     71      INTEGER                              , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4) 
     72      INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4) 
    7673      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
    7774      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
     
    7976      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
    8077      ! 
    81       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
    82       INTEGER  ::   ik   
    83       REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar 
    84       REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      - 
    85       REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      - 
    86       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz 
    87       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 
    88       !!---------------------------------------------------------------------- 
    89       ! 
    90       IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd') 
    91       ! 
    92       CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz ) 
     78      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices   
     79      REAL(wp) ::   z2dtt, ztra                              ! local scalar 
     80      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      - 
     81      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      - 
     82      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 
     83      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz 
     84      !!---------------------------------------------------------------------- 
     85      ! 
     86      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct') 
     87      ! 
     88      CALL wrk_alloc( jpi,jpj,jpk,   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 
    9389      ! 
    9490      IF( kt == kit000 )  THEN 
    9591         IF(lwp) WRITE(numout,*) 
    96          IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype 
     92         IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype 
    9793         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    98          ! 
    99          l_trd = .FALSE. 
    100          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    10194      ENDIF 
     95      ! 
     96      l_trd = .FALSE. 
     97      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    10298      ! 
    10399      IF( l_trd )  THEN 
    104100         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
    105          ztrdx(:,:,:) = 0.e0   ;    ztrdy(:,:,:) = 0.e0   ;   ztrdz(:,:,:) = 0.e0 
     101         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp 
    106102      ENDIF 
    107103      ! 
    108       zwi(:,:,:) = 0.e0 ;  
    109       ! 
     104      !                                         ! surface & bottom value : flux set to zero one for all 
     105      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp                ! except at the surface in linear free surface case 
     106      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp 
     107      ! 
     108      zwi(:,:,:) = 0._wp         
    110109      !                                                          ! =========== 
    111110      DO jn = 1, kjpt                                            ! tracer loop 
    112111         !                                                       ! =========== 
    113          ! 1. Bottom and k=1 value : flux set to zero 
    114          ! ---------------------------------- 
    115          zwx(:,:,jpk) = 0.e0    ;    zwz(:,:,jpk) = 0.e0 
    116          zwy(:,:,jpk) = 0.e0    ;    zwi(:,:,jpk) = 0.e0 
    117            
    118          zwz(:,:,1  ) = 0._wp 
    119          ! 2. upstream advection with initial mass fluxes & intermediate update 
    120          ! -------------------------------------------------------------------- 
    121          ! upstream tracer flux in the i and j direction 
     112         ! 
     113         !        !==  upstream advection with initial mass fluxes & intermediate update  ==! 
     114         !                    !* upstream tracer flux in the i and j direction  
    122115         DO jk = 1, jpkm1 
    123116            DO jj = 1, jpjm1 
     
    133126            END DO 
    134127         END DO 
    135  
    136          ! upstream tracer flux in the k direction 
    137          ! Interior value 
    138          DO jk = 2, jpkm1 
     128         !                    !* upstream tracer flux in the k direction *! 
     129         DO jk = 2, jpkm1         ! Interior value ( multiplied by wmask) 
    139130            DO jj = 1, jpj 
    140131               DO ji = 1, jpi 
     
    145136            END DO 
    146137         END DO 
    147          ! Surface value 
    148          IF( lk_vvl ) THEN    
    149             IF ( ln_isfcav ) THEN 
    150                DO jj = 1, jpj 
    151                   DO ji = 1, jpi 
    152                      zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable 
    153                   END DO 
    154                END DO 
    155             ELSE 
    156                zwz(:,:,1) = 0.e0          ! volume variable 
    157             END IF 
    158          ELSE                 
    159             IF ( ln_isfcav ) THEN 
     138         !                     
     139         IF(.NOT.lk_vvl ) THEN   ! top ocean value (only in linear free surface as zwz has been w-masked) 
     140            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface 
    160141               DO jj = 1, jpj 
    161142                  DO ji = 1, jpi 
     
    163144                  END DO 
    164145               END DO    
    165             ELSE 
    166                zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
    167             END IF 
     146            ELSE                             ! no cavities: only at the ocean surface 
     147               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     148            ENDIF 
    168149         ENDIF 
    169  
    170          ! total advective trend 
    171          DO jk = 1, jpkm1 
     150         !                
     151         DO jk = 1, jpkm1     !* trend and after field with monotonic scheme 
    172152            z2dtt = p2dt(jk) 
    173153            DO jj = 2, jpjm1 
    174154               DO ji = fs_2, fs_jpim1   ! vector opt. 
    175                   zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    176155                  ! total intermediate advective trends 
    177                   ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    178                      &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    179                      &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
     156                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     157                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     158                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    180159                  ! update and guess with monotonic sheme 
     160!!gm why tmask added in the two following lines ???    the mask is done in tranxt ! 
    181161                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra   * tmask(ji,jj,jk) 
    182162                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     
    184164            END DO 
    185165         END DO 
    186          !                             ! Lateral boundary conditions on zwi  (unchanged sign) 
    187          CALL lbc_lnk( zwi, 'T', 1. )   
    188  
    189          !                                 ! trend diagnostics (contribution of upstream fluxes) 
    190          IF( l_trd )  THEN  
    191             ! store intermediate advective trends 
     166         CALL lbc_lnk( zwi, 'T', 1. )  ! Lateral boundary conditions on zwi  (unchanged sign) 
     167         !                 
     168         IF( l_trd )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
    192169            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:) 
    193170         END IF 
    194          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     171         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    195172         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    196173           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    197174           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    198175         ENDIF 
    199  
    200          ! 3. antidiffusive flux : high order minus low order 
    201          ! -------------------------------------------------- 
    202          ! antidiffusive flux on i and j 
    203          DO jk = 1, jpkm1 
    204             DO jj = 1, jpjm1 
    205                DO ji = 1, fs_jpim1   ! vector opt. 
    206                   zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk) 
    207                   zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk) 
    208                END DO 
    209             END DO 
    210          END DO 
    211        
    212          ! antidiffusive flux on k 
    213          ! Interior value 
    214          DO jk = 2, jpkm1                     
    215             DO jj = 1, jpj 
    216                DO ji = 1, jpi 
    217                   zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 
    218                END DO 
    219             END DO 
    220          END DO 
    221          ! surface value 
    222          IF ( ln_isfcav ) THEN 
    223             DO jj = 1, jpj 
    224                DO ji = 1, jpi 
    225                   zwz(ji,jj,mikt(ji,jj)) = 0.e0 
    226                END DO 
    227             END DO 
    228          ELSE 
    229             zwz(:,:,1) = 0.e0 
    230          END IF 
     176         ! 
     177         ! 
     178         !        !==  anti-diffusive flux : high order minus low order  ==! 
     179         ! 
     180         SELECT CASE( kn_fct_h )         !* horizontal anti-diffusive fluxes 
     181         ! 
     182         CASE(  2  )                         ! 2nd order centered 
     183            DO jk = 1, jpkm1 
     184               DO jj = 1, jpjm1 
     185                  DO ji = 1, fs_jpim1   ! vector opt. 
     186                     zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk) 
     187                     zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk) 
     188                  END DO 
     189               END DO 
     190            END DO 
     191            ! 
     192         CASE(  4  )                         ! 4th order centered 
     193            zltu(:,:,jpk) = 0._wp                            ! Bottom value : flux set to zero 
     194            zltv(:,:,jpk) = 0._wp 
     195            DO jk = 1, jpkm1                                 ! Laplacian 
     196               DO jj = 1, jpjm1                                   ! First derivative (gradient) 
     197                  DO ji = 1, fs_jpim1   ! vector opt. 
     198                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     199                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     200                  END DO 
     201               END DO 
     202               DO jj = 2, jpjm1                                   !  
     203                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     204                     zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6 
     205                     zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6 
     206                  END DO 
     207               END DO 
     208            END DO 
     209            ! 
     210            CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
     211            ! 
     212            DO jk = 1, jpkm1                                 ! Horizontal advective fluxes 
     213               DO jj = 1, jpjm1 
     214                  DO ji = 1, fs_jpim1   ! vector opt. 
     215                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points 
     216                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) 
     217                     !                                                  ! C4 minus upstream advective fluxes  
     218                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 
     219                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 
     220                  END DO 
     221               END DO 
     222            END DO          
     223            ! 
     224         CASE(  41 )                         ! 4th order centered       ==>>   !!gm coding attempt   need to be tested 
     225            ztu(:,:,jpk) = 0._wp                             ! Bottom value : flux set to zero 
     226            ztv(:,:,jpk) = 0._wp 
     227            DO jk = 1, jpkm1                                 ! gradient 
     228               DO jj = 1, jpjm1                                   ! First derivative (gradient) 
     229                  DO ji = 1, fs_jpim1   ! vector opt. 
     230                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     231                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     232                  END DO 
     233               END DO 
     234            END DO 
     235            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn) 
     236            ! 
     237            DO jk = 1, jpkm1                                 ! Horizontal advective fluxes 
     238               DO jj = 2, jpjm1 
     239                  DO ji = 2, fs_jpim1   ! vector opt. 
     240                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points (x2) 
     241                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) 
     242                     !                                                  ! C4 interpolation of T at u- & v-points (x2) 
     243                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) ) 
     244                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) ) 
     245                     !                                                  ! C4 minus upstream advective fluxes  
     246                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 
     247                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 
     248                  END DO 
     249               END DO 
     250            END DO 
     251            ! 
     252         END SELECT 
     253         !                                !* vertical anti-diffusive fluxes 
     254         SELECT CASE( kn_fct_v )                ! Interior values (w-masked) 
     255         ! 
     256         CASE(  2  )                                  ! 2nd order centered 
     257            DO jk = 2, jpkm1     
     258               DO jj = 2, jpjm1 
     259                  DO ji = fs_2, fs_jpim1 
     260                     zwz(ji,jj,jk) =  ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
     261                                       - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
     262                  END DO 
     263               END DO 
     264            END DO 
     265            ! 
     266         CASE(  4  )                                  ! 4th order COMPACT 
     267            !     
     268            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! COMPACT interpolation of T at w-point 
     269            ! 
     270            DO jk = 2, jpkm1 
     271               DO jj = 2, jpjm1 
     272                  DO ji = fs_2, fs_jpim1 
     273                     zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
     274                  END DO 
     275               END DO 
     276            END DO 
     277            ! 
     278         END SELECT 
     279         !                                      ! top ocean value: high order = upstream  ==>>  zwz=0 
     280         zwz(:,:, 1 ) = 0._wp                   ! only ocean surface as interior zwz values have been w-masked 
     281         ! 
    231282         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    232283         CALL lbc_lnk( zwz, 'W',  1. ) 
    233284 
    234          ! 4. monotonicity algorithm 
    235          ! ------------------------- 
     285         !        !==  monotonicity algorithm  ==! 
     286         ! 
    236287         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 
    237288 
    238289 
    239          ! 5. final trend with corrected fluxes 
    240          ! ------------------------------------ 
     290         !        !==  final trend with corrected fluxes  ==! 
     291         ! 
    241292         DO jk = 1, jpkm1 
    242293            DO jj = 2, jpjm1 
    243294               DO ji = fs_2, fs_jpim1   ! vector opt.   
    244                   zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    245                   ! total advective trends 
    246                   ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    247                      &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    248                      &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
    249                   ! add them to the general tracer trends 
    250                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 
    251                END DO 
    252             END DO 
    253          END DO 
    254  
    255          !                                 ! trend diagnostics (contribution of upstream fluxes) 
    256          IF( l_trd )  THEN  
     295                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     296                     &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     297                     &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) & 
     298                     &                                / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     299               END DO 
     300            END DO 
     301         END DO 
     302         ! 
     303         IF( l_trd ) THEN                 ! trend diagnostics (contribution of upstream fluxes) 
    257304            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed 
    258305            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
    259306            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
    260              
    261             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )    
    262             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )   
    263             CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )  
     307            ! 
     308            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 
     309            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
     310            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
     311            ! 
     312            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz ) 
    264313         END IF 
    265314         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     
    271320      END DO 
    272321      ! 
    273                    CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz ) 
    274       IF( l_trd )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
    275       ! 
    276       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd') 
    277       ! 
    278    END SUBROUTINE tra_adv_tvd 
    279  
    280  
    281    SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    282       &                                       ptb, ptn, pta, kjpt ) 
    283       !!---------------------------------------------------------------------- 
    284       !!                  ***  ROUTINE tra_adv_tvd_zts  *** 
     322      CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 
     323      ! 
     324      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct') 
     325      ! 
     326   END SUBROUTINE tra_adv_fct 
     327 
     328 
     329   SUBROUTINE tra_adv_fct_zts( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
     330      &                                                  ptb, ptn, pta, kjpt, kn_fct_zts ) 
     331      !!---------------------------------------------------------------------- 
     332      !!                  ***  ROUTINE tra_adv_fct_zts  *** 
    285333      !!  
    286334      !! **  Purpose :   Compute the now trend due to total advection of  
     
    296344      !!             - save the trends  
    297345      !!---------------------------------------------------------------------- 
    298       USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace 
    299       ! 
    300346      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    301347      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    302348      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    303349      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
     350      INTEGER                              , INTENT(in   ) ::   kn_fct_zts      ! number of number of vertical sub-timesteps 
    304351      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
    305352      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
     
    310357      REAL(wp), DIMENSION( jpk )                           ::   zr_p2dt         ! reciprocal of tracer timestep 
    311358      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices   
    312       INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection 
    313359      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps 
    314360      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps 
    315361      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection 
    316       REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar 
     362      REAL(wp) ::   z2dtt, ztra              ! local scalar 
    317363      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      - 
    318364      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      - 
    319       REAL(wp), POINTER, DIMENSION(:,:  ) :: zwx_sav , zwy_sav 
    320       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts 
    321       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 
    322       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs 
    323       !!---------------------------------------------------------------------- 
    324       ! 
    325       IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd_zts') 
    326       ! 
    327       CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav ) 
    328       CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts ) 
    329       CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs ) 
     365      REAL(wp), POINTER, DIMENSION(:,:  )   ::  zwx_sav , zwy_sav 
     366      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav 
     367      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdx, ztrdy, ztrdz 
     368      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrs 
     369      !!---------------------------------------------------------------------- 
     370      ! 
     371      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct_zts') 
     372      ! 
     373      CALL wrk_alloc( jpi,jpj,            zwx_sav, zwy_sav ) 
     374      CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
     375      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1, ztrs ) 
    330376      ! 
    331377      IF( kt == kit000 )  THEN 
    332378         IF(lwp) WRITE(numout,*) 
    333          IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype 
     379         IF(lwp) WRITE(numout,*) 'tra_adv_fct_zts : 2nd order FCT scheme with ', kn_fct_zts, ' vertical sub-timestep on ', cdtype 
    334380         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    335381      ENDIF 
    336382      ! 
    337383      l_trd = .FALSE. 
    338       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
     384      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    339385      ! 
    340386      IF( l_trd )  THEN 
    341          CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
     387         CALL wrk_alloc( jpi,jpj,jpk,  ztrdx, ztrdy, ztrdz ) 
    342388         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp 
    343389      ENDIF 
    344390      ! 
    345391      zwi(:,:,:) = 0._wp 
    346       z_rzts = 1._wp / REAL( jnzts, wp ) 
     392      z_rzts = 1._wp / REAL( kn_fct_zts, wp ) 
    347393      zr_p2dt(:) = 1._wp / p2dt(:) 
    348394      ! 
     
    373419 
    374420         ! upstream tracer flux in the k direction 
    375          ! Interior value 
    376          DO jk = 2, jpkm1 
     421         DO jk = 2, jpkm1         ! Interior value 
    377422            DO jj = 1, jpj 
    378423               DO ji = 1, jpi 
    379424                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    380425                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    381                   zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 
    382                END DO 
    383             END DO 
    384          END DO 
    385          ! Surface value 
    386          IF( lk_vvl ) THEN 
    387             IF ( ln_isfcav ) THEN 
     426                  zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
     427               END DO 
     428            END DO 
     429         END DO 
     430         !                       ! top value 
     431         IF( lk_vvl ) THEN             ! variable volume: only k=1 as zwz is multiplied by wmask 
     432            zwz(:,:, 1 ) = 0._wp 
     433         ELSE                          ! linear free surface 
     434            IF( ln_isfcav ) THEN             ! ice-shelf cavities 
    388435               DO jj = 1, jpj 
    389436                  DO ji = 1, jpi 
    390                      zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable +    isf 
    391                   END DO 
    392                END DO 
    393             ELSE 
    394                zwz(:,:,1) = 0.e0                              ! volume variable + no isf 
    395             END IF 
    396          ELSE 
    397             IF ( ln_isfcav ) THEN 
    398                DO jj = 1, jpj 
    399                   DO ji = 1, jpi 
    400                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface +    isf 
    401                   END DO 
    402                END DO 
    403             ELSE 
    404                zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf 
    405             END IF 
     437                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)  
     438                  END DO 
     439               END DO    
     440            ELSE                             ! standard case 
     441               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     442            ENDIF 
    406443         ENDIF 
    407  
    408          ! total advective trend 
    409          DO jk = 1, jpkm1 
     444         ! 
     445         DO jk = 1, jpkm1         ! total advective trend 
    410446            z2dtt = p2dt(jk) 
    411447            DO jj = 2, jpjm1 
    412448               DO ji = fs_2, fs_jpim1   ! vector opt. 
    413                   zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    414449                  ! total intermediate advective trends 
    415                   ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    416                      &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    417                      &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
     450                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     451                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     452                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    418453                  ! update and guess with monotonic sheme 
    419454                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra 
     
    422457            END DO 
    423458         END DO 
    424          !                             ! Lateral boundary conditions on zwi  (unchanged sign) 
    425          CALL lbc_lnk( zwi, 'T', 1. )   
    426  
    427          !                                 ! trend diagnostics (contribution of upstream fluxes) 
    428          IF( l_trd )  THEN  
    429             ! store intermediate advective trends 
     459         !                            
     460         CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign) 
     461         !                 
     462         IF( l_trd )  THEN                ! trend diagnostics (contribution of upstream fluxes) 
    430463            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:) 
    431464         END IF 
    432          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     465         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    433466         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    434467           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
     
    436469         ENDIF 
    437470 
    438          ! 3. antidiffusive flux : high order minus low order 
    439          ! -------------------------------------------------- 
    440          ! antidiffusive flux on i and j 
    441  
    442  
    443          DO jk = 1, jpkm1 
    444  
     471         ! 3. anti-diffusive flux : high order minus low order 
     472         ! --------------------------------------------------- 
     473 
     474         DO jk = 1, jpkm1                    !* horizontal anti-diffusive fluxes 
     475            ! 
    445476            DO jj = 1, jpjm1 
    446477               DO ji = 1, fs_jpim1   ! vector opt. 
    447478                  zwx_sav(ji,jj) = zwx(ji,jj,jk) 
    448479                  zwy_sav(ji,jj) = zwy(ji,jj,jk) 
    449  
     480                  ! 
    450481                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) 
    451482                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) 
    452483               END DO 
    453484            END DO 
    454  
    455             DO jj = 2, jpjm1         ! partial horizontal divergence 
     485            ! 
     486            DO jj = 2, jpjm1                    ! partial horizontal divergence 
    456487               DO ji = fs_2, fs_jpim1 
    457488                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
     
    459490               END DO 
    460491            END DO 
    461  
     492            ! 
    462493            DO jj = 1, jpjm1 
    463494               DO ji = 1, fs_jpim1   ! vector opt. 
    464                   zwx(ji,jj,jk) = zwx(ji,jj,jk)  - zwx_sav(ji,jj) 
    465                   zwy(ji,jj,jk) = zwy(ji,jj,jk)  - zwy_sav(ji,jj) 
     495                  zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj) 
     496                  zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj) 
    466497               END DO 
    467498            END DO 
    468499         END DO 
    469500       
    470          ! antidiffusive flux on k 
    471          zwz(:,:,1) = 0._wp        ! Surface value 
    472          zwz_sav(:,:,:) = zwz(:,:,:) 
    473          ! 
    474          ztrs(:,:,:,1) = ptb(:,:,:,jn) 
    475          zwzts(:,:,:) = 0._wp 
    476  
    477          DO jl = 1, jnzts                   ! Start of sub timestepping loop 
    478  
    479             IF( jl == 1 ) THEN              ! Euler forward to kick things off 
    480               jtb = 1   ;   jtn = 1   ;   jta = 2 
    481               zts(:) = p2dt(:) * z_rzts 
    482               jtaken = MOD( jnzts + 1 , 2)  ! Toggle to collect every second flux 
    483                                             ! starting at jl =1 if jnzts is odd;  
    484                                             ! starting at jl =2 otherwise 
    485             ELSEIF( jl == 2 ) THEN          ! First leapfrog step 
    486               jtb = 1   ;   jtn = 2   ;   jta = 3 
    487               zts(:) = 2._wp * p2dt(:) * z_rzts 
    488             ELSE                            ! Shuffle pointers for subsequent leapfrog steps 
    489               jtb = MOD(jtb,3) + 1 
    490               jtn = MOD(jtn,3) + 1 
    491               jta = MOD(jta,3) + 1 
     501         !                                !* vertical anti-diffusive flux 
     502         zwz_sav(:,:,:)   = zwz(:,:,:) 
     503         ztrs   (:,:,:,1) = ptb(:,:,:,jn) 
     504         zwzts  (:,:,:)   = 0._wp 
     505         IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp    ! surface value set to zero in vvl case 
     506         ! 
     507         DO jl = 1, kn_fct_zts                  ! Start of sub timestepping loop 
     508            ! 
     509            IF( jl == 1 ) THEN                        ! Euler forward to kick things off 
     510               jtb = 1   ;   jtn = 1   ;   jta = 2 
     511               zts(:) = p2dt(:) * z_rzts 
     512               jtaken = MOD( kn_fct_zts + 1 , 2)            ! Toggle to collect every second flux 
     513               !                                            ! starting at jl =1 if kn_fct_zts is odd;  
     514               !                                            ! starting at jl =2 otherwise 
     515            ELSEIF( jl == 2 ) THEN                    ! First leapfrog step 
     516               jtb = 1   ;   jtn = 2   ;   jta = 3 
     517               zts(:) = 2._wp * p2dt(:) * z_rzts 
     518            ELSE                                      ! Shuffle pointers for subsequent leapfrog steps 
     519               jtb = MOD(jtb,3) + 1 
     520               jtn = MOD(jtn,3) + 1 
     521               jta = MOD(jta,3) + 1 
    492522            ENDIF 
    493             DO jk = 2, jpkm1          ! Interior value 
     523            DO jk = 2, jpkm1                          ! interior value 
    494524               DO jj = 2, jpjm1 
    495525                  DO ji = fs_2, fs_jpim1 
    496                      zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) 
    497                      IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk)           ! Accumulate time-weighted vertcal flux 
    498                   END DO 
    499                END DO 
    500             END DO 
    501  
     526                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) * wmask(ji,jj,jk) 
     527                     IF( jtaken == 0 )   zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk) * zts(jk)    ! Accumulate time-weighted vertcal flux 
     528                  END DO 
     529               END DO 
     530            END DO 
     531            IF(.NOT.lk_vvl ) THEN                    ! top value (only in linear free surface case) 
     532               IF( ln_isfcav ) THEN                      ! ice-shelf cavities 
     533                  DO jj = 1, jpj 
     534                     DO ji = 1, jpi 
     535                        zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     536                     END DO 
     537                  END DO    
     538               ELSE                                      ! standard case 
     539                  zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     540               ENDIF 
     541            ENDIF 
     542            ! 
    502543            jtaken = MOD( jtaken + 1 , 2 ) 
    503  
    504             DO jk = 2, jpkm1          ! Interior value 
     544            ! 
     545            DO jk = 2, jpkm1                             ! total advective trends 
    505546               DO jj = 2, jpjm1 
    506547                  DO ji = fs_2, fs_jpim1 
    507                      zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    508                      ! total advective trends 
    509                      ztra = - zbtr * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
    510                      ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra 
    511                   END DO 
    512                END DO 
    513             END DO 
    514  
     548                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb)                                                 & 
     549                        &               - zts(jk) * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
     550                        &                         / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     551                  END DO 
     552               END DO 
     553            END DO 
     554            ! 
    515555         END DO 
    516556 
     
    518558            DO jj = 2, jpjm1 
    519559               DO ji = fs_2, fs_jpim1 
    520                   zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) 
     560                  zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk) 
    521561               END DO 
    522562            END DO 
     
    535575            DO jj = 2, jpjm1 
    536576               DO ji = fs_2, fs_jpim1   ! vector opt.   
    537                   zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    538                   ! total advective trends 
    539                   ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    540                      &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    541                      &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
    542                   ! add them to the general tracer trends 
    543                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     577                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )       & 
     578                     &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   )   & 
     579                     &                                / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    544580               END DO 
    545581            END DO 
     
    551587            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
    552588            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
    553              
     589            ! 
    554590            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )    
    555591            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )   
    556592            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )  
     593            ! 
     594            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz ) 
    557595         END IF 
    558596         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     
    564602      END DO 
    565603      ! 
    566                    CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts ) 
    567                    CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs ) 
    568                    CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav ) 
    569       IF( l_trd )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
    570       ! 
    571       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd_zts') 
    572       ! 
    573    END SUBROUTINE tra_adv_tvd_zts 
     604      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav ) 
     605      CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
     606      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs ) 
     607      ! 
     608      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts') 
     609      ! 
     610   END SUBROUTINE tra_adv_fct_zts 
    574611 
    575612 
     
    683720   END SUBROUTINE nonosc 
    684721 
     722 
     723   SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 
     724      !!---------------------------------------------------------------------- 
     725      !!                  ***  ROUTINE interp_4th_cpt  *** 
     726      !!  
     727      !! **  Purpose :   Compute the interpolation of tracer at w-point 
     728      !! 
     729      !! **  Method  :   4th order compact interpolation 
     730      !!---------------------------------------------------------------------- 
     731      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! now tracer fields 
     732      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! now tracer field interpolated at w-pts 
     733      ! 
     734      INTEGER :: ji, jj, jk   ! dummy loop integers 
     735      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 
     736      !!---------------------------------------------------------------------- 
     737       
     738      DO jk = 3, jpkm1        !==  build the three diagonal matrix  ==! 
     739         DO jj = 1, jpj 
     740            DO ji = 1, jpi 
     741               zwd (ji,jj,jk) = 4._wp 
     742               zwi (ji,jj,jk) = 1._wp 
     743               zws (ji,jj,jk) = 1._wp 
     744               zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
     745               ! 
     746               IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom 
     747                  zwd (ji,jj,jk) = 1._wp 
     748                  zwi (ji,jj,jk) = 0._wp 
     749                  zws (ji,jj,jk) = 0._wp 
     750                  zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )     
     751               ENDIF 
     752            END DO 
     753         END DO 
     754      END DO 
     755      ! 
     756      jk=2                                            ! Switch to second order centered at top 
     757      DO jj=1,jpj 
     758         DO ji=1,jpi 
     759            zwd (ji,jj,jk) = 1._wp 
     760            zwi (ji,jj,jk) = 0._wp 
     761            zws (ji,jj,jk) = 0._wp 
     762            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
     763         END DO 
     764      END DO    
     765      ! 
     766      !                       !==  tridiagonal solve  ==! 
     767      DO jj = 1, jpj                ! first recurrence 
     768         DO ji = 1, jpi 
     769            zwt(ji,jj,2) = zwd(ji,jj,2) 
     770         END DO 
     771      END DO 
     772      DO jk = 3, jpkm1 
     773         DO jj = 1, jpj 
     774            DO ji = 1, jpi 
     775               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
     776            END DO 
     777         END DO 
     778      END DO 
     779      ! 
     780      DO jj = 1, jpj                ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
     781         DO ji = 1, jpi 
     782            pt_out(ji,jj,2) = zwrm(ji,jj,2) 
     783         END DO 
     784      END DO 
     785      DO jk = 3, jpkm1 
     786         DO jj = 1, jpj 
     787            DO ji = 1, jpi 
     788               pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
     789            END DO 
     790         END DO 
     791      END DO 
     792 
     793      DO jj = 1, jpj                ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 
     794         DO ji = 1, jpi 
     795            pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
     796         END DO 
     797      END DO 
     798      DO jk = jpk-2, 2, -1 
     799         DO jj = 1, jpj 
     800            DO ji = 1, jpi 
     801               pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
     802            END DO 
     803         END DO 
     804      END DO 
     805      !     
     806   END SUBROUTINE interp_4th_cpt 
     807    
    685808   !!====================================================================== 
    686 END MODULE traadv_tvd 
     809END MODULE traadv_fct 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90

    r5737 r5770  
    1 MODULE traadv_muscl 
     1MODULE traadv_mus 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  traadv_muscl  *** 
     3   !!                       ***  MODULE  traadv_mus  *** 
    44   !! Ocean  tracers:  horizontal & vertical advective trend 
    55   !!====================================================================== 
     
    99   !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport 
    1010   !!            3.4  !  2012-06  (P. Oddo, M. Vichi) include the upstream where needed 
    11    !!---------------------------------------------------------------------- 
    12  
    13    !!---------------------------------------------------------------------- 
    14    !!   tra_adv_muscl : update the tracer trend with the horizontal 
     11   !!            3.7  !  2015-09  (G. Madec) add the ice-shelf cavities boundary condition 
     12   !!---------------------------------------------------------------------- 
     13 
     14   !!---------------------------------------------------------------------- 
     15   !!   tra_adv_mus   : update the tracer trend with the horizontal 
    1516   !!                   and vertical advection trends using MUSCL scheme 
    1617   !!---------------------------------------------------------------------- 
     
    3435   PRIVATE 
    3536 
    36    PUBLIC   tra_adv_muscl   ! routine called by traadv.F90 
     37   PUBLIC   tra_adv_mus   ! routine called by traadv.F90 
    3738    
    3839   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits 
     
    5051CONTAINS 
    5152 
    52    SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn,   & 
    53       &                                                ptb, pta, kjpt, ld_msc_ups ) 
     53   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn,             & 
     54      &                                              ptb, pta, kjpt, ld_msc_ups ) 
    5455      !!---------------------------------------------------------------------- 
    55       !!                    ***  ROUTINE tra_adv_muscl  *** 
    56       !! 
    57       !! ** Purpose :   Compute the now trend due to total advection of T and  
    58       !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for 
    59       !!      Conservation Laws) and add it to the general tracer trend. 
     56      !!                    ***  ROUTINE tra_adv_mus  *** 
     57      !! 
     58      !! ** Purpose :   Compute the now trend due to total advection of tracers 
     59      !!              using a MUSCL scheme (Monotone Upstream-centered Scheme for 
     60      !!              Conservation Laws) and add it to the general tracer trend. 
    6061      !! 
    6162      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries 
     63      !!              ld_msc_ups=T :  
    6264      !! 
    6365      !! ** Action  : - update (ta,sa) with the now advective tracer trends 
    64       !!              - save trends  
     66      !!              - send trends to trdtra module for further diagnostcs 
    6567      !! 
    6668      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation 
     
    7779      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
    7880      ! 
    79       INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices 
    80       INTEGER  ::   ierr                      ! local integer 
    81       REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars 
    82       REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      - 
    83       REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      - 
     81      INTEGER  ::   ji, jj, jk, jn       ! dummy loop indices 
     82      INTEGER  ::   ierr                 ! local integer 
     83      REAL(wp) ::   zu, z0u, zzwx, zw    ! local scalars 
     84      REAL(wp) ::   zv, z0v, zzwy, z0w   !   -      - 
     85      REAL(wp) ::   zdt, zalpha          !   -      - 
    8486      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zslpx, zslpy   ! 3D workspace 
    8587      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx  , zwy     ! -      -  
    8688      !!---------------------------------------------------------------------- 
    8789      ! 
    88       IF( nn_timing == 1 )  CALL timing_start('tra_adv_muscl') 
    89       ! 
    90       CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy ) 
     90      IF( nn_timing == 1 )  CALL timing_start('tra_adv_mus') 
     91      ! 
     92      CALL wrk_alloc( jpi,jpj,jpk,  zslpx, zslpy, zwx, zwy ) 
    9193      ! 
    9294      IF( kt == kit000 )  THEN 
     
    9799         IF(lwp) WRITE(numout,*) 
    98100         ! 
    99          ! 
    100          IF( ld_msc_ups ) THEN 
    101             IF( .NOT. ALLOCATED( upsmsk ) )  THEN 
    102                 ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 
    103                 IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate upsmsk array') 
    104             ENDIF 
     101         ! Upstream / MUSCL scheme indicator 
     102         ! 
     103         ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr ) 
     104         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed 
     105         ! 
     106         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked) 
     107            ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 
    105108            upsmsk(:,:) = 0._wp                             ! not upstream by default 
    106          ENDIF 
    107  
    108          IF( .NOT. ALLOCATED( xind ) ) THEN 
    109              ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr ) 
    110              IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate zind array') 
    111          ENDIF 
    112          ! 
    113  
    114          ! 
    115          ! Upstream / MUSCL scheme indicator 
    116          ! ------------------------------------ 
    117 !!gm  useless 
    118          xind(:,:,:) = 1._wp                             ! set equal to 1 where up-stream is not needed 
    119 !!gm 
    120          ! 
    121          IF( ld_msc_ups )  THEN 
     109            ! 
    122110            DO jk = 1, jpkm1 
    123111               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed 
    124112                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows) 
    125                   &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 near some straits 
     113                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area 
    126114            END DO 
    127115         ENDIF  
     
    172160               END DO 
    173161           END DO 
    174          END DO             ! interior values 
    175  
     162         END DO 
     163         ! 
    176164         !                                             !-- MUSCL horizontal advective fluxes 
    177165         DO jk = 1, jpkm1                                     ! interior values 
     
    196184            END DO 
    197185         END DO 
    198          !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign) 
    199          CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. ) 
    200          ! 
    201          ! Tracer flux divergence at t-point added to the general trend 
    202          DO jk = 1, jpkm1 
     186         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign) 
     187         ! 
     188         DO jk = 1, jpkm1         ! Tracer flux divergence at t-point added to the general trend 
    203189            DO jj = 2, jpjm1       
    204190               DO ji = fs_2, fs_jpim1   ! vector opt. 
    205                   zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    206                   ! horizontal advective trends 
    207                   ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    208                   &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
    209                   ! add it to the general tracer trends 
    210                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     191                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       & 
     192                  &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
     193                  &                                   / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    211194               END DO 
    212195           END DO 
     
    226209         ! II. Vertical advective fluxes 
    227210         ! ----------------------------- 
    228          !                                             !-- first guess of the slopes 
    229          zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0    ! surface & bottom boundary conditions 
     211         !                                !-- first guess of the slopes 
     212         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions 
     213         zwx(:,:,jpk) = 0._wp                   ! surface & bottom boundary conditions 
    230214         DO jk = 2, jpkm1                                     ! interior values 
    231215            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 
    232216         END DO 
    233217 
    234          !                                             !-- Slopes of tracer 
    235          zslpx(:,:,1) = 0.e0                                  ! surface values 
    236          DO jk = 2, jpkm1                                     ! interior value 
     218         !                                !-- Slopes of tracer 
     219         zslpx(:,:,1) = 0._wp                   ! surface values 
     220         DO jk = 2, jpkm1                       ! interior value 
    237221            DO jj = 1, jpj 
    238222               DO ji = 1, jpi 
     
    242226            END DO 
    243227         END DO 
    244          !                                             !-- Slopes limitation 
    245          DO jk = 2, jpkm1                                     ! interior values 
     228         !                                !-- Slopes limitation 
     229         DO jk = 2, jpkm1                       ! interior values 
    246230            DO jj = 1, jpj 
    247231               DO ji = 1, jpi 
     
    252236            END DO 
    253237         END DO 
    254          !                                             !-- vertical advective flux 
    255          !                                                    ! surface values  (bottom already set to zero) 
    256          IF( lk_vvl ) THEN    ;   zwx(:,:, 1 ) = 0.e0                      !  variable volume 
    257          ELSE                 ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
    258          ENDIF  
    259          ! 
    260          DO jk = 1, jpkm1                                     ! interior values 
     238         !                                !-- vertical advective flux 
     239         DO jk = 1, jpkm1                       ! interior values 
    261240            zdt  = p2dt(jk) 
    262241            DO jj = 2, jpjm1       
    263242               DO ji = fs_2, fs_jpim1   ! vector opt. 
    264                   zbtr = 1. / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 
    265243                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
    266244                  zalpha = 0.5 + z0w 
    267                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr  
     245                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt  / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 
    268246                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
    269247                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
    270                   zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     248                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
    271249               END DO  
    272250            END DO 
    273251         END DO 
    274  
     252         !                                      ! top values  (bottom already set to zero) 
     253         IF( lk_vvl ) THEN                            !  variable volume 
     254            zwx(:,:, 1 ) = 0._wp                            ! k=1 only as zwx has been multiplied by wmask 
     255         ELSE                                         ! linear free surface  
     256            IF( ln_isfcav ) THEN                            ! ice-shelf cavities (top of the ocean) 
     257               DO jj = 1, jpj 
     258                  DO ji = 1, jpi 
     259                     zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 
     260                  END DO 
     261               END DO    
     262            ELSE                                             ! no cavities: only at the ocean surface 
     263               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     264            ENDIF 
     265         ENDIF 
     266         ! 
    275267         DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend 
    276268            DO jj = 2, jpjm1       
    277269               DO ji = fs_2, fs_jpim1   ! vector opt. 
    278                   zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    279                   ! vertical advective trends  
    280                   ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
    281                   ! add it to the general tracer trends 
    282                   pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) + ztra 
     270                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    283271               END DO 
    284272            END DO 
     
    291279      END DO 
    292280      ! 
    293       CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy ) 
    294       ! 
    295       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl') 
    296       ! 
    297    END SUBROUTINE tra_adv_muscl 
     281      CALL wrk_dealloc( jpi,jpj,jpk,  zslpx, zslpy, zwx, zwy ) 
     282      ! 
     283      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_mus') 
     284      ! 
     285   END SUBROUTINE tra_adv_mus 
    298286 
    299287   !!====================================================================== 
    300 END MODULE traadv_muscl 
     288END MODULE traadv_mus 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r5737 r5770  
    102102         IF(lwp) WRITE(numout,*) 
    103103      ENDIF 
     104      ! 
    104105      l_trd = .FALSE. 
    105106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
     
    130131      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
    131132      !! 
    132       INTEGER  :: ji, jj, jk, jn   ! dummy loop indices 
    133       REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars 
    134       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd 
     133      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     134      REAL(wp) ::   ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars 
     135      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zfu, zfc, zfd 
    135136      !---------------------------------------------------------------------- 
    136137      ! 
     
    139140      DO jn = 1, kjpt                                            ! tracer loop 
    140141         !                                                       ! =========== 
    141          zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0   
    142          zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0      
    143          !                                                   
    144          DO jk = 1, jpkm1                                 
    145             !                                              
    146             !--- Computation of the ustream and downstream value of the tracer and the mask 
     142         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp  
     143         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp    
     144         ! 
     145!!gm why not using a SHIFT instruction... 
     146         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask 
    147147            DO jj = 2, jpjm1 
    148148               DO ji = fs_2, fs_jpim1   ! vector opt. 
    149                   ! Upstream in the x-direction for the tracer 
    150                   zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) 
    151                   ! Downstream in the x-direction for the tracer 
    152                   zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) 
     149                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer 
     150                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer 
    153151               END DO 
    154152            END DO 
     
    159157         ! Horizontal advective fluxes 
    160158         ! --------------------------- 
    161          ! 
    162159         DO jk = 1, jpkm1                              
    163160            DO jj = 2, jpjm1 
     
    380377      ! 
    381378      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    382       REAL(wp) ::   zbtr , ztra      ! local scalars 
    383379      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz 
    384380      !!---------------------------------------------------------------------- 
    385381      ! 
    386       CALL wrk_alloc( jpi, jpj, jpk, zwz ) 
     382      CALL wrk_alloc( jpi,jpj,jpk,   zwz ) 
     383      ! 
     384      !                          ! surface & bottom values  
     385      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all 
     386                     zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface 
     387      ! 
    387388      !                                                          ! =========== 
    388389      DO jn = 1, kjpt                                            ! tracer loop 
    389390         !                                                       ! =========== 
    390          ! 1. Bottom value : flux set to zero 
    391          zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero 
    392          ! 
    393          !                                 ! Surface value 
    394          IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero 
    395          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface 
     391         ! 
     392         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux) 
     393            DO jj = 2, jpjm1 
     394               DO ji = fs_2, fs_jpim1   ! vector opt. 
     395                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 
     396               END DO 
     397            END DO 
     398         END DO 
     399         IF(.NOT.lk_vvl ) THEN               !* top value   (only in linear free surf. as zwz is multiplied by wmask) 
     400            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean) 
     401               DO jj = 1, jpj 
     402                  DO ji = 1, jpi 
     403                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     404                  END DO 
     405               END DO    
     406            ELSE                                   ! no ice-shelf cavities (only ocean surface) 
     407               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 
     408            ENDIF 
    396409         ENDIF 
    397410         ! 
    398          DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point 
     411         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==! 
    399412            DO jj = 2, jpjm1 
    400413               DO ji = fs_2, fs_jpim1   ! vector opt. 
    401                   zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) 
    402                END DO 
    403             END DO 
    404          END DO 
    405          ! 
    406          DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==! 
    407             DO jj = 2, jpjm1 
    408                DO ji = fs_2, fs_jpim1   ! vector opt. 
    409                   zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    410                   ! k- vertical advective trends  
    411                   ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )  
    412                   ! added to the general tracer trends 
    413                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     414                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
     415                     &                                / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    414416               END DO 
    415417            END DO 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r5737 r5770  
    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  
    1819   USE trdtra         ! trends manager: tracers  
    1920   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
     
    3839#  include "vectopt_loop_substitute.h90" 
    3940   !!---------------------------------------------------------------------- 
    40    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     41   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    4142   !! $Id$ 
    4243   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    4445CONTAINS 
    4546 
    46    SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    47       &                                       ptb, ptn, pta, kjpt ) 
     47   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn,          & 
     48      &                                                ptb, ptn, pta, kjpt, kn_ubs_v ) 
    4849      !!---------------------------------------------------------------------- 
    4950      !!                  ***  ROUTINE tra_adv_ubs  *** 
     
    5253      !!      and add it to the general trend of passive tracer equations. 
    5354      !! 
    54       !! ** Method  :   The upstream biased scheme (UBS) is based on a 3rd order 
     55      !! ** Method  :   The 3rd order Upstream Biased Scheme (UBS) is based on an 
    5556      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 
    5657      !!      It is only used in the horizontal direction. 
     
    6162      !!      where zltu is the second derivative of the before temperature field: 
    6263      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 
    63       !!      This results in a dissipatively dominant (i.e. hyper-diffusive)  
     64      !!        This results in a dissipatively dominant (i.e. hyper-diffusive)  
    6465      !!      truncation error. The overall performance of the advection scheme  
    6566      !!      is similar to that reported in (Farrow and Stevens, 1995).  
    66       !!      For stability reasons, the first term of the fluxes which corresponds 
     67      !!        For stability reasons, the first term of the fluxes which corresponds 
    6768      !!      to a second order centered scheme is evaluated using the now velocity  
    6869      !!      (centered in time) while the second term which is the diffusive part  
    6970      !!      of the scheme, is evaluated using the before velocity (forward in time).  
    7071      !!      Note that UBS is not positive. Do not use it on passive tracers. 
    71       !!                On the vertical, the advection is evaluated using a TVD scheme, 
    72       !!      as the UBS have been found to be too diffusive. 
     72      !!                On the vertical, the advection is evaluated using a FCT scheme, 
     73      !!      as the UBS have been found to be too diffusive.  
     74!!gm  !!                kn_ubs_v argument (not coded for the moment) 
     75      !!      controles whether the FCT is based on a 2nd order centrered scheme (kn_ubs_v=2)  
     76      !!      or on a 4th order compact scheme (kn_ubs_v=4). 
    7377      !! 
    7478      !! ** Action : - update (pta) with the now advective tracer trends 
     
    8185      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    8286      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
     87      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers 
    8388      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
    8489      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components 
     
    95100      IF( nn_timing == 1 )  CALL timing_start('tra_adv_ubs') 
    96101      ! 
    97       CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 
     102      CALL wrk_alloc( jpi,jpj,jpk,  ztu, ztv, zltu, zltv, zti, ztw ) 
    98103      ! 
    99104      IF( kt == kit000 )  THEN 
     
    106111      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    107112      ! 
     113      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp     ! Bottom value : set to zero one for all 
     114      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp 
     115      IF( lk_vvl )   ztw(:,:, 1 ) = 0._wp                   ! surface value: set to zero only in vvl case 
     116      ! 
    108117      !                                                          ! =========== 
    109118      DO jn = 1, kjpt                                            ! tracer loop 
    110119         !                                                       ! =========== 
    111          ! 1. Bottom value : flux set to zero 
    112          ! ---------------------------------- 
    113          zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0 
    114120         !                                               
    115          DO jk = 1, jpkm1                                 ! Horizontal slab 
    116             !                                    
    117             !  Laplacian 
    118             DO jj = 1, jpjm1            ! First derivative (gradient) 
     121         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==! 
     122            DO jj = 1, jpjm1              ! First derivative (masked gradient) 
    119123               DO ji = 1, fs_jpim1   ! vector opt. 
    120124                  zeeu = e2_e1u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk) 
     
    124128               END DO 
    125129            END DO 
    126             DO jj = 2, jpjm1            ! Second derivative (divergence) 
     130            DO jj = 2, jpjm1              ! Second derivative (divergence) 
    127131               DO ji = fs_2, fs_jpim1   ! vector opt. 
    128                   zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 
     132                  zcoef = 1._wp / ( 6._wp * fse3t(ji,jj,jk) ) 
    129133                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef 
    130134                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
     
    132136            END DO 
    133137            !                                     
    134          END DO                                           ! End of slab          
     138         END DO          
    135139         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    136  
    137140         !     
    138          !  Horizontal advective fluxes                
    139          DO jk = 1, jpkm1                                 ! Horizontal slab 
     141         DO jk = 1, jpkm1        !==  Horizontal advective fluxes  ==!     (UBS) 
    140142            DO jj = 1, jpjm1 
    141143               DO ji = 1, fs_jpim1   ! vector opt. 
    142                   ! upstream transport (x2) 
    143                   zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
     144                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! upstream transport (x2) 
    144145                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
    145146                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
    146147                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
    147                   ! 2nd order centered advective fluxes (x2) 
     148                  !                                                  ! 2nd order centered advective fluxes (x2) 
    148149                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) ) 
    149150                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) ) 
    150                   ! UBS advective fluxes 
     151                  !                                                  ! UBS advective fluxes 
    151152                  ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 
    152153                  ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 
    153154               END DO 
    154155            END DO 
    155          END DO                                           ! End of slab          
    156  
    157          zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends 
    158  
    159          DO jk = 1, jpkm1                 ! Horizontal advective trends 
     156         END DO          
     157         ! 
     158         zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update 
     159         ! 
     160         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==! 
    160161            DO jj = 2, jpjm1 
    161162               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    166167            END DO 
    167168            !                                              
    168          END DO                                           !   End of slab 
    169  
    170          ! Horizontal trend used in tra_adv_ztvd subroutine 
    171          zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) 
    172  
     169         END DO 
     170         ! 
     171         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case 
     172         !                                            ! and/or in trend diagnostic (l_trd=T)  
    173173         !                 
    174174         IF( l_trd ) THEN                  ! trend diagnostics 
     
    181181            IF( jn == jp_sal )  str_adv(:) = ptr_sj( ztv(:,:,:) ) 
    182182         ENDIF 
    183           
    184          ! TVD scheme for the vertical direction   
    185          ! ---------------------- 
    186          IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag. 
    187  
    188          !  Bottom value : flux set to zero 
    189          ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0 
    190  
    191          ! Surface value 
    192          IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero 
    193          ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface  
    194          ENDIF 
    195          !  upstream advection with initial mass fluxes & intermediate update 
    196          ! ------------------------------------------------------------------- 
    197          ! Interior value 
    198          DO jk = 2, jpkm1 
    199             DO jj = 1, jpj 
    200                DO ji = 1, jpi 
    201                    zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    202                    zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    203                    ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) 
    204                END DO 
    205             END DO 
    206          END DO  
    207          ! update and guess with monotonic sheme 
    208          DO jk = 1, jpkm1 
    209             z2dtt = p2dt(jk) 
    210             DO jj = 2, jpjm1 
    211                DO ji = fs_2, fs_jpim1   ! vector opt. 
    212                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    213                   ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 
    214                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak  
    215                   zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
    216                END DO 
    217             END DO 
    218          END DO 
    219          ! 
    220          CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
    221  
    222          !  antidiffusive flux : high order minus low order 
    223          ztw(:,:,1) = 0.e0       ! Surface value 
    224          DO jk = 2, jpkm1        ! Interior value 
    225             DO jj = 1, jpj 
    226                DO ji = 1, jpi 
    227                   ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk) 
    228                END DO 
    229             END DO 
    230          END DO 
    231          ! 
    232          CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm 
    233  
    234          !  final trend with corrected fluxes 
    235          DO jk = 1, jpkm1 
     183         ! 
     184         !                       !== vertical advective trend  ==! 
     185         ! 
     186         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme 
     187         ! 
     188         CASE(  2  )                   ! 2nd order FCT  
     189            !          
     190            IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag. 
     191            ! 
     192            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==! 
     193            DO jk = 2, jpkm1                 ! Interior value (w-masked) 
     194               DO jj = 1, jpj 
     195                  DO ji = 1, jpi 
     196                     zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     197                     zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     198                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk) 
     199                  END DO 
     200               END DO 
     201            END DO  
     202            IF(.NOT.lk_vvl ) THEN            ! top ocean value (only in linear free surface as ztw has been w-masked) 
     203               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface 
     204                  DO jj = 1, jpj 
     205                     DO ji = 1, jpi 
     206                        ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     207                     END DO 
     208                  END DO    
     209               ELSE                                ! no cavities: only at the ocean surface 
     210                  ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     211               ENDIF 
     212            ENDIF 
     213            ! 
     214            DO jk = 1, jpkm1           !* trend and after field with monotonic scheme 
     215               z2dtt = p2dt(jk) 
     216               DO jj = 2, jpjm1 
     217                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     218                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
     219                     pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak  
     220                     zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
     221                  END DO 
     222               END DO 
     223            END DO 
     224            CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
     225            ! 
     226            !                          !*  anti-diffusive flux : high order minus low order 
     227            DO jk = 2, jpkm1        ! Interior value  (w-masked) 
     228               DO jj = 1, jpj 
     229                  DO ji = 1, jpi 
     230                     ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
     231                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk) 
     232                  END DO 
     233               END DO 
     234            END DO 
     235            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0 
     236            IF(.NOT.lk_vvl )   ztw(:,:, 1 ) = 0._wp      ! only ocean surface as interior zwz values have been w-masked 
     237            ! 
     238            CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm 
     239            ! 
     240         CASE(  4  )                               ! 4th order COMPACT 
     241            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point 
     242            DO jk = 2, jpkm1 
     243               DO jj = 2, jpjm1 
     244                  DO ji = fs_2, fs_jpim1 
     245                     ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
     246                  END DO 
     247               END DO 
     248            END DO 
     249            IF(.NOT.lk_vvl )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work 
     250            ! 
     251         END SELECT 
     252         ! 
     253         DO jk = 1, jpkm1        !  final trend with corrected fluxes 
    236254            DO jj = 2, jpjm1  
    237255               DO ji = fs_2, fs_jpim1   ! vector opt.    
    238                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    239                   ! k- vertical advective trends   
    240                   ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) 
    241                   ! added to the general tracer trends 
    242                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    243                END DO 
    244             END DO 
    245          END DO 
    246  
    247          !  Save the final vertical advective trends 
    248          IF( l_trd )  THEN                        ! vertical advective trend diagnostics 
     256                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
     257               END DO 
     258            END DO 
     259         END DO 
     260         ! 
     261         IF( l_trd )  THEN       ! vertical advective trend diagnostics 
    249262            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 
    250263               DO jj = 2, jpjm1 
    251264                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    252                      zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    253                      z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr 
    254                      zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn 
     265                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          & 
     266                        &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   & 
     267                        &                              / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    255268                  END DO 
    256269               END DO 
     
    261274      END DO 
    262275      ! 
    263       CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 
     276      CALL wrk_dealloc( jpi,jpj,jpk,  ztu, ztv, zltu, zltv, zti, ztw ) 
    264277      ! 
    265278      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs') 
     
    294307      IF( nn_timing == 1 )  CALL timing_start('nonosc_z') 
    295308      ! 
    296       CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo ) 
     309      CALL wrk_alloc( jpi,jpj,jpk,  zbetup, zbetdo ) 
    297310      ! 
    298311      zbig  = 1.e+40_wp 
    299312      zrtrn = 1.e-15_wp 
    300313      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp 
    301  
     314      ! 
    302315      ! Search local extrema 
    303316      ! -------------------- 
    304       ! large negative value (-zbig) inside land 
     317      !                    ! large negative value (-zbig) inside land 
    305318      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 
    306319      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 
    307       ! search maximum in neighbourhood 
    308       DO jk = 1, jpkm1 
     320      ! 
     321      DO jk = 1, jpkm1     ! search maximum in neighbourhood 
    309322         ikm1 = MAX(jk-1,1) 
    310323         DO jj = 2, jpjm1 
     
    316329         END DO 
    317330      END DO 
    318       ! large positive value (+zbig) inside land 
     331      !                    ! large positive value (+zbig) inside land 
    319332      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 
    320333      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 
    321       ! search minimum in neighbourhood 
    322       DO jk = 1, jpkm1 
     334      ! 
     335      DO jk = 1, jpkm1     ! search minimum in neighbourhood 
    323336         ikm1 = MAX(jk-1,1) 
    324337         DO jj = 2, jpjm1 
     
    330343         END DO 
    331344      END DO 
    332  
    333       ! restore masked values to zero 
     345      !                    ! restore masked values to zero 
    334346      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 
    335347      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 
    336  
    337  
    338       ! 2. Positive and negative part of fluxes and beta terms 
    339       ! ------------------------------------------------------ 
    340  
     348      ! 
     349      ! Positive and negative part of fluxes and beta terms 
     350      ! --------------------------------------------------- 
    341351      DO jk = 1, jpkm1 
    342352         z2dtt = p2dt(jk) 
     
    347357               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
    348358               ! up & down beta terms 
    349                zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 
     359               zbt = e1e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 
    350360               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 
    351361               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt 
     
    353363         END DO 
    354364      END DO 
     365      ! 
    355366      ! monotonic flux in the k direction, i.e. pcc 
    356367      ! ------------------------------------------- 
     
    366377      END DO 
    367378      ! 
    368       CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo ) 
     379      CALL wrk_dealloc( jpi,jpj,jpk,  zbetup, zbetdo ) 
    369380      ! 
    370381      IF( nn_timing == 1 )  CALL timing_stop('nonosc_z') 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r5766 r5770  
    7373      SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend 
    7474      ! 
    75       CASE ( n_lap   )                                   ! laplacian: iso-level operator 
     75      CASE ( np_lap   )                                  ! laplacian: iso-level operator 
    7676         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   ) 
    7777         ! 
    78       CASE ( n_lap_i )                                   ! laplacian: standard iso-neutral operator (Madec) 
     78      CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec) 
    7979         CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
    8080         ! 
    81       CASE ( n_lap_it )                                  ! laplacian: triad iso-neutral operator (griffies) 
     81      CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies) 
    8282         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
    8383         ! 
    84       CASE ( n_blp , n_blp_i , n_blp_it )                ! bilaplacian: iso-level & iso-neutral operators 
     84      CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators 
    8585         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf ) 
    8686      END SELECT 
     
    126126      IF( ln_traldf_blp )   ioptio = ioptio + 1 
    127127      IF( ioptio >  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on tracer' ) 
    128       IF( ioptio == 0   )   nldf = n_no_ldf     ! No lateral diffusion 
     128      IF( ioptio == 0   )   nldf = np_no_ldf     ! No lateral diffusion 
    129129      ioptio = 0 
    130130      IF( ln_traldf_lev )   ioptio = ioptio + 1 
     
    137137      IF( ln_traldf_lap ) THEN         ! laplacian operator 
    138138         IF ( ln_zco ) THEN               ! z-coordinate 
    139             IF ( ln_traldf_lev   )   nldf = n_lap     ! iso-level = horizontal (no rotation) 
    140             IF ( ln_traldf_hor   )   nldf = n_lap     ! iso-level = horizontal (no rotation) 
    141             IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard  (   rotation) 
    142             IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad     (   rotation) 
     139            IF ( ln_traldf_lev   )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     140            IF ( ln_traldf_hor   )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     141            IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation) 
     142            IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation) 
    143143         ENDIF 
    144144         IF ( ln_zps ) THEN               ! z-coordinate with partial step 
    145             IF ( ln_traldf_lev   )   ierr = 1         ! iso-level not allowed  
    146             IF ( ln_traldf_hor   )   nldf = n_lap     ! horizontal (no rotation) 
    147             IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation) 
    148             IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation) 
     145            IF ( ln_traldf_lev   )   ierr = 1          ! iso-level not allowed  
     146            IF ( ln_traldf_hor   )   nldf = np_lap     ! horizontal            (no rotation) 
     147            IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard    (rotation) 
     148            IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad        (rotation) 
    149149         ENDIF 
    150150         IF ( ln_sco ) THEN               ! s-coordinate 
    151             IF ( ln_traldf_lev   )   nldf = n_lap     ! iso-level  (no rotation) 
    152             IF ( ln_traldf_hor   )   nldf = n_lap_i   ! horizontal (   rotation)    
    153             IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation) 
    154             IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation) 
     151            IF ( ln_traldf_lev   )   nldf = np_lap     ! iso-level              (no rotation) 
     152            IF ( ln_traldf_hor   )   nldf = np_lap_i   ! horizontal             (   rotation) 
     153            IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation) 
     154            IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation) 
    155155         ENDIF 
    156156      ENDIF 
     
    158158      IF( ln_traldf_blp ) THEN         ! bilaplacian operator 
    159159         IF ( ln_zco ) THEN               ! z-coordinate 
    160             IF ( ln_traldf_lev   )   nldf = n_blp     ! iso-level = horizontal (no rotation) 
    161             IF ( ln_traldf_hor   )   nldf = n_blp     ! iso-level = horizontal (no rotation) 
    162             IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
    163             IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     160            IF ( ln_traldf_lev   )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
     161            IF ( ln_traldf_hor   )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
     162            IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
     163            IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation) 
    164164         ENDIF 
    165165         IF ( ln_zps ) THEN               ! z-coordinate with partial step 
    166             IF ( ln_traldf_lev   )   ierr = 1         ! iso-level not allowed  
    167             IF ( ln_traldf_hor   )   nldf = n_blp     ! horizontal (no rotation) 
    168             IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
    169             IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     166            IF ( ln_traldf_lev   )   ierr = 1          ! iso-level not allowed  
     167            IF ( ln_traldf_hor   )   nldf = np_blp     ! horizontal            (no rotation) 
     168            IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
     169            IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation) 
    170170         ENDIF 
    171171         IF ( ln_sco ) THEN               ! s-coordinate 
    172             IF ( ln_traldf_lev   )   nldf = n_blp     ! iso-level  (no rotation) 
    173             IF ( ln_traldf_hor   )   nldf = n_blp_it  ! horizontal (   rotation)       !!gm   a checker.... 
    174             IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
    175             IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     172            IF ( ln_traldf_lev   )   nldf = np_blp     ! iso-level              (no rotation) 
     173            IF ( ln_traldf_hor   )   nldf = np_blp_it  ! horizontal             (   rotation) 
     174            IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
     175            IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation) 
    176176         ENDIF 
    177177      ENDIF 
     
    181181           &            CALL ctl_stop( '          eddy induced velocity on tracers requires isopycnal',    & 
    182182           &                                                                    ' laplacian diffusion' ) 
    183       IF(  nldf == n_lap_i .OR. nldf == n_lap_it .OR. & 
    184          & nldf == n_blp_i .OR. nldf == n_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required  
     183      IF(  nldf == np_lap_i .OR. nldf == np_lap_it .OR. & 
     184         & nldf == np_blp_i .OR. nldf == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required  
    185185      ! 
    186186      IF(lwp) THEN 
    187187         WRITE(numout,*) 
    188          IF( nldf == n_no_ldf )   WRITE(numout,*) '          NO lateral diffusion' 
    189          IF( nldf == n_lap    )   WRITE(numout,*) '          laplacian iso-level operator' 
    190          IF( nldf == n_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)' 
    191          IF( nldf == n_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)' 
    192          IF( nldf == n_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator' 
    193          IF( nldf == n_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)' 
    194          IF( nldf == n_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)' 
     188         IF( nldf == np_no_ldf )   WRITE(numout,*) '          NO lateral diffusion' 
     189         IF( nldf == np_lap    )   WRITE(numout,*) '          laplacian iso-level operator' 
     190         IF( nldf == np_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)' 
     191         IF( nldf == np_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)' 
     192         IF( nldf == np_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator' 
     193         IF( nldf == np_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)' 
     194         IF( nldf == np_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)' 
    195195      ENDIF 
    196196      ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_blp.F90

    r5758 r5770  
    88 
    99   !!---------------------------------------------------------------------- 
    10    !!   traldf_blp   : update the tracer trend with the bilaplacian lateral mixing trend 
     10   !!   traldf_blp    : update the tracer trend with the bilaplacian lateral mixing trend 
    1111   !!---------------------------------------------------------------------- 
    12    USE oce           ! ocean dynamics and active tracers 
    13    USE dom_oce       ! ocean space and time domain 
    14    USE phycst        ! physical constants 
    15    USE trc_oce       ! share passive tracers/Ocean variables 
    16    USE zdf_oce       ! ocean vertical physics 
    17    USE ldftra        ! lateral physics: eddy diffusivity 
    18    USE ldfslp        ! lateral physics: iso-neutral slopes 
    19    USE traldf_lap    ! lateral diffusion (Standard operator)         (tra_ldf_lap   routine) 
    20    USE traldf_iso    ! lateral diffusion (Standard operator)         (tra_ldf_iso   routine) 
    21    USE traldf_triad  ! lateral diffusion (Standard operator)         (tra_ldf_triad routine) 
    22    USE diaptr        ! poleward transport diagnostics 
    23    USE zpshde        ! partial step: hor. derivative     (zps_hde routine) 
     12   USE oce            ! ocean dynamics and active tracers 
     13   USE dom_oce        ! ocean space and time domain 
     14   USE phycst         ! physical constants 
     15   USE trc_oce        ! share passive tracers/Ocean variables 
     16   USE zdf_oce        ! ocean vertical physics 
     17   USE ldftra         ! lateral physics: eddy diffusivity 
     18   USE ldfslp         ! lateral physics: iso-neutral slopes 
     19   USE traldf_lap     ! lateral diffusion (Standard operator)         (tra_ldf_lap   routine) 
     20   USE traldf_iso     ! lateral diffusion (Standard operator)         (tra_ldf_iso   routine) 
     21   USE traldf_triad   ! lateral diffusion (Standard operator)         (tra_ldf_triad routine) 
     22   USE diaptr         ! poleward transport diagnostics 
     23   USE zpshde         ! partial step: hor. derivative     (zps_hde routine) 
    2424   ! 
    2525   USE in_out_manager ! I/O manager 
     
    3333   PRIVATE 
    3434 
    35    PUBLIC   tra_ldf_blp         ! routine called by traldf.F90 
     35   PUBLIC   tra_ldf_blp   ! routine called by traldf.F90 
    3636 
    37    !                                    ! Flag to control the type of lateral diffusive operator 
    38    INTEGER, PARAMETER, PUBLIC ::   n_no_ldf = 00                      ! without operator (i.e. no lateral diffusive trend) 
    39    !               !!      laplacian    !    bilaplacian   ! 
    40    INTEGER, PARAMETER, PUBLIC ::   n_lap    = 10   ,   n_blp    = 20  ! iso-level operator 
    41    INTEGER, PARAMETER, PUBLIC ::   n_lap_i  = 11   ,   n_blp_i  = 21  ! standard iso-neutral or geopotential operator  
    42    INTEGER, PARAMETER, PUBLIC ::   n_lap_it = 12   ,   n_blp_it = 22  ! triad    iso-neutral or geopotential operator 
     37   !                      ! Flag to control the type of lateral diffusive operator 
     38   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral diffusive trend) 
     39   !                          !!      laplacian     !    bilaplacian    ! 
     40   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator 
     41   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator 
     42   INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator 
    4343 
    4444   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels 
     
    9494         WRITE(numout,*) 
    9595         SELECT CASE ( kldf ) 
    96          CASE ( n_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype 
    97          CASE ( n_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)' 
    98          CASE ( n_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)' 
     96         CASE ( np_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype 
     97         CASE ( np_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)' 
     98         CASE ( np_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)' 
    9999         END SELECT 
    100100         WRITE(numout,*) '~~~~~~~~~~~' 
     
    105105      SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==! 
    106106      ! 
    107       CASE ( n_blp    )                ! iso-level bilaplacian 
     107      CASE ( np_blp    )               ! iso-level bilaplacian 
    108108         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 ) 
    109109         ! 
    110       CASE ( n_blp_i  )                ! rotated   bilaplacian : standard operator (Madec) 
     110      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec) 
    111111         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 
    112112         ! 
    113       CASE ( n_blp_it )                ! rotated  bilaplacian : triad operator (griffies) 
     113      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies) 
    114114         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 
    115115         ! 
     
    126126      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==! 
    127127      ! 
    128       CASE ( n_blp    )                ! iso-level bilaplacian 
     128      CASE ( np_blp    )               ! iso-level bilaplacian 
    129129         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 ) 
    130130         ! 
    131       CASE ( n_blp_i  )                ! rotated   bilaplacian : standard operator (Madec) 
     131      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec) 
    132132         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 
    133133         ! 
    134       CASE ( n_blp_it )                ! rotated  bilaplacian : triad operator (griffies) 
     134      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies) 
    135135         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 
    136136         ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r5758 r5770  
    8383      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme  
    8484      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra,            tsb, tsa, jpts )  !   implicit scheme  
    85       CASE ( -1 )                                       ! esopa: test all possibility with control print 
    86          CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts ) 
    87          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf0 - Ta: ', mask1=tmask,               & 
    88          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    89          CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra,            tsb, tsa, jpts )  
    90          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf1 - Ta: ', mask1=tmask,               & 
    91          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    9285      END SELECT 
    9386!!gm WHY here !   and I don't like that ! 
     
    149142            &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' ) 
    150143 
    151       ! Test: esopa 
    152       IF( lk_esopa )    nzdf = -1                      ! All schemes used 
    153  
    154144      IF(lwp) THEN 
    155145         WRITE(numout,*) 
    156146         WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme' 
    157147         WRITE(numout,*) '~~~~~~~~~~~' 
    158          IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used' 
    159148         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme' 
    160149         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme' 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90

    r5758 r5770  
    800800      END IF 
    801801 
    802       IF( nn_cla == 1 )   CALL ctl_warn( '      You set n_cla = 1. Note that the Mixed-Layer diagnostics  ',   & 
    803          &                               '      are not exact along the corresponding straits.            ') 
    804  
    805802      !                                   ! allocate trdmxl arrays 
    806803      IF( trd_mxl_alloc()    /= 0 )   CALL ctl_stop( 'STOP', 'trd_mxl_init : unable to allocate trdmxl     arrays' ) 
     
    809806 
    810807 
    811  
    812        nkstp     = nit000 - 1              ! current time step indicator initialization 
     808      nkstp     = nit000 - 1              ! current time step indicator initialization 
    813809 
    814810 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r5758 r5770  
    3333   !!             -   ! 2013-06  (I. Epicoco, S. Mocavero, CMCC) nemo_northcomms: setup avoiding MPI communication  
    3434   !!            3.7  ! 2014-12  (G. Madec) suppression of cross land advection option 
    35    !!             -   ! 2014-12  (G. Madec) remove KPP scheme 
     35   !!             -   ! 2014-12  (G. Madec) remove KPP scheme and cross-land advection (cla) 
    3636   !!---------------------------------------------------------------------- 
    3737 
     
    4646   !!---------------------------------------------------------------------- 
    4747   USE step_oce        ! module used in the ocean time stepping module 
    48    USE cla             ! cross land advection               (tra_cla routine) 
    4948   USE domcfg          ! domain configuration               (dom_cfg routine) 
    5049   USE mppini          ! shared/distributed memory setting (mpp_init routine) 
     
    458457 
    459458      !                                      ! Misc. options 
    460       IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_init       ! Cross Land Advection 
    461459                            CALL sto_par_init   ! Stochastic parametrization 
    462460      IF( ln_sto_eos     )  CALL sto_pts_init   ! RRandom T/S fluctuations 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/step.F90

    r5758 r5770  
    109109      IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
    110110                         CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice) 
    111 CALL FLUSH    ( numout ) 
     111 
    112112      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    113113      ! Update stochastic parameters and random T/S fluctuations 
    114114      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    115115                         CALL sto_par( kstp )          ! Stochastic parameters 
    116 CALL FLUSH    ( numout ) 
    117116 
    118117      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    124123                         CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 
    125124                         CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency 
    126 CALL FLUSH    ( numout ) 
     125 
    127126      ! 
    128127      !  VERTICAL PHYSICS 
     
    137136         avmv(:,:,:) = rn_avm0 * wvmask(:,:,:) 
    138137      ENDIF 
    139 CALL FLUSH    ( numout ) 
     138 
    140139      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths 
    141140         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO 
     
    152151      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' ) 
    153152      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' ) 
    154 CALL FLUSH    ( numout ) 
    155153      ! 
    156154      !  LATERAL  PHYSICS 
     
    179177      !                                               ! eddy diffusivity coeff. and/or eiv coeff. 
    180178      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp )  
    181 write(*,*) 'after ldf_slp' 
    182 CALL FLUSH    ( numout ) 
    183179 
    184180      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    188184      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    189185                         CALL wzv           ( kstp )  ! now cross-level velocity  
    190 write(*,*) 'after wzv' 
    191 CALL FLUSH    ( numout ) 
    192186 
    193187      IF( lk_dynspg_ts ) THEN  
     
    231225                                  CALL wzv           ( kstp )  ! now cross-level velocity  
    232226      ENDIF 
    233 write(*,*) 'after wzv 2' 
    234 CALL FLUSH    ( numout ) 
    235227 
    236228      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    253245                         CALL trc_stp( kstp )         ! time-stepping 
    254246#endif 
    255 write(*,*) 'end dyn ' 
    256 CALL FLUSH    ( numout ) 
    257247 
    258248      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    270260      IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends 
    271261                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection 
    272 write(*,*) 'before tra_ldf' 
    273 CALL FLUSH    ( numout ) 
    274262                             CALL tra_ldf    ( kstp )       ! lateral mixing 
    275 write(*,*) 'after tra_ldf' 
    276 CALL FLUSH    ( numout ) 
    277263 
    278264!!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 
     
    316302                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
    317303      ENDIF 
    318 write(*,*) 'after tra_nxt' 
    319 CALL FLUSH    ( numout ) 
    320304 
    321305      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    360344      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
    361345      ! 
    362 write(*,*) 'after dom_vvl' 
    363 CALL FLUSH    ( numout ) 
    364  
    365346 
    366347!!gm : This does not only concern the dynamics ==>>> add a new title 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90

    r5766 r5770  
    44   !! Ocean passive tracers:  advection trend  
    55   !!============================================================================== 
    6    !! History :  2.0  !  05-11  (G. Madec)  Original code 
    7    !!            3.0  !  10-06  (C. Ethe)   Adapted to passive tracers 
     6   !! History :  2.0  !  2005-11  (G. Madec)  Original code 
     7   !!            3.0  !  2010-06  (C. Ethe)   Adapted to passive tracers 
     8   !!            3.7  !  2014-05  (G. Madec, C. Ethe)  Add 2nd/4th order cases for CEN and FCT schemes  
    89   !!---------------------------------------------------------------------- 
    910#if defined key_top 
     
    1112   !!   'key_top'                                                TOP models 
    1213   !!---------------------------------------------------------------------- 
    13    !!   trc_adv      : compute ocean tracer advection trend 
    14    !!   trc_adv_ini  : control the different options of advection scheme 
    15    !!---------------------------------------------------------------------- 
    16    USE oce_trc         ! ocean dynamics and active tracers 
    17    USE trc             ! ocean passive tracers variables 
    18    USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine) 
    19    USE traadv_tvd      ! TVD      scheme           (tra_adv_tvd    routine) 
    20    USE traadv_muscl    ! MUSCL    scheme           (tra_adv_muscl  routine) 
    21    USE traadv_muscl2   ! MUSCL2   scheme           (tra_adv_muscl2 routine) 
    22    USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine) 
    23    USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine) 
    24    USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine) 
    25    USE ldftra          ! lateral diffusion coefficient on tracers 
    26    USE prtctl_trc      ! Print control 
     14   !!   trc_adv       : compute ocean tracer advection trend 
     15   !!   trc_adv_ini   : control the different options of advection scheme 
     16   !!---------------------------------------------------------------------- 
     17   USE oce_trc        ! ocean dynamics and active tracers 
     18   USE trc            ! ocean passive tracers variables 
     19   USE traadv_cen     ! centered scheme           (tra_adv_cen  routine) 
     20   USE traadv_fct     ! FCT      scheme           (tra_adv_fct  routine) 
     21   USE traadv_mus     ! MUSCL    scheme           (tra_adv_mus  routine) 
     22   USE traadv_ubs     ! UBS      scheme           (tra_adv_ubs  routine) 
     23   USE traadv_qck     ! QUICKEST scheme           (tra_adv_qck  routine) 
     24   USE traadv_mle     ! ML eddy induced velocity  (tra_adv_mle  routine) 
     25   USE ldftra         ! lateral diffusion coefficient on tracers 
     26   USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces 
     27   ! 
     28   USE prtctl_trc     ! Print control 
    2729 
    2830   IMPLICIT NONE 
     
    3335   PUBLIC   trc_adv_ini   
    3436 
    35    !                                        !!: ** Advection (namtrc_adv) ** 
    36    LOGICAL , PUBLIC ::   ln_trcadv_cen2      ! 2nd order centered scheme flag 
    37    LOGICAL , PUBLIC ::   ln_trcadv_tvd       ! TVD scheme flag 
    38    LOGICAL , PUBLIC ::   ln_trcadv_muscl     ! MUSCL scheme flag 
    39    LOGICAL , PUBLIC ::   ln_trcadv_muscl2    ! MUSCL2 scheme flag 
    40    LOGICAL , PUBLIC ::   ln_trcadv_ubs       ! UBS scheme flag 
    41    LOGICAL , PUBLIC ::   ln_trcadv_qck       ! QUICKEST scheme flag 
    42    LOGICAL , PUBLIC ::   ln_trcadv_msc_ups   ! use upstream scheme within muscl 
    43  
    44    INTEGER ::   nadv   ! choice of the type of advection scheme 
     37   !                            !!* Namelist namtrc_adv * 
     38   LOGICAL ::   ln_trcadv_cen    ! centered scheme flag 
     39   INTEGER ::      nn_cen_h, nn_cen_v   ! =2/4 : horizontal and vertical choices of the order of CEN scheme 
     40   LOGICAL ::   ln_trcadv_fct    ! FCT scheme flag 
     41   INTEGER ::      nn_fct_h, nn_fct_v   ! =2/4 : horizontal and vertical choices of the order of FCT scheme 
     42   INTEGER ::      nn_fct_zts           ! >=1 : 2nd order FCT with vertical sub-timestepping 
     43   LOGICAL ::   ln_trcadv_mus    ! MUSCL scheme flag 
     44   LOGICAL ::      ln_mus_ups           ! use upstream scheme in vivcinity of river mouths 
     45   LOGICAL ::   ln_trcadv_ubs    ! UBS scheme flag 
     46   INTEGER ::      nn_ubs_v             ! =2/4 : vertical choice of the order of UBS scheme 
     47   LOGICAL ::   ln_trcadv_qck    ! QUICKEST scheme flag 
     48 
     49   !                                        ! choices of advection scheme: 
     50   INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection 
     51   INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme 
     52   INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme 
     53   INTEGER, PARAMETER ::   np_FCT_zts = 3   ! 2nd order FCT scheme with vertical sub-timestepping 
     54   INTEGER, PARAMETER ::   np_MUS     = 4   ! MUSCL scheme 
     55   INTEGER, PARAMETER ::   np_UBS     = 5   ! 3rd order Upstream Biased Scheme 
     56   INTEGER, PARAMETER ::   np_QCK     = 6   ! QUICK scheme 
     57 
     58   INTEGER ::              nadv             ! chosen advection scheme 
     59   ! 
    4560   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dt  ! vertical profile time-step, = 2 rdttra 
    4661   !                                                    ! except at nitrrc000 (=rdttra) if neuler=0 
     
    5065#  include "vectopt_loop_substitute.h90" 
    5166   !!---------------------------------------------------------------------- 
    52    !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     67   !! NEMO/TOP 3.7 , NEMO Consortium (2015) 
    5368   !! $Id$  
    5469   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     
    6075      !!                  ***  ROUTINE trc_adv_alloc  *** 
    6176      !!---------------------------------------------------------------------- 
    62  
     77      ! 
    6378      ALLOCATE( r2dt(jpk), STAT=trc_adv_alloc ) 
    64  
     79      ! 
    6580      IF( trc_adv_alloc /= 0 ) CALL ctl_warn('trc_adv_alloc : failed to allocate array.') 
    66  
     81      ! 
    6782   END FUNCTION trc_adv_alloc 
    6883 
     
    7893      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    7994      ! 
    80       INTEGER ::   jk  
     95      INTEGER ::   jk   ! dummy loop index 
    8196      CHARACTER (len=22) ::   charout 
    8297      REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn  ! effective velocity 
     
    93108      ENDIF 
    94109      !                                               !==  effective transport  ==! 
    95       zun(:,:,jpk) = 0._wp                                                       ! no transport trough the bottom 
    96       zvn(:,:,jpk) = 0._wp 
    97       zwn(:,:,jpk) = 0._wp 
    98110      DO jk = 1, jpkm1 
    99111         zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                   ! eulerian transport 
     
    112124      IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )  ! add the mle transport 
    113125      ! 
    114       ! 
    115       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==! 
    116       ! 
    117       CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, nittrc000, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered 
    118       CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD  
    119       CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra, ln_trcadv_msc_ups )   !  MUSCL  
    120       CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2  
    121       CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS  
    122       CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  QUICKEST  
     126      zun(:,:,jpk) = 0._wp                                                       ! no transport trough the bottom 
     127      zvn(:,:,jpk) = 0._wp 
     128      zwn(:,:,jpk) = 0._wp 
     129      ! 
     130      ! 
     131      SELECT CASE ( nadv )                      !==  compute advection trend and add it to general trend  ==! 
     132      ! 
     133      CASE ( np_CEN )                                    ! Centered : 2nd / 4th order 
     134         CALL tra_adv_cen    ( kt, nittrc000,'TRC',       zun, zvn, zwn     , trn, tra, jptra, nn_cen_h, nn_cen_v ) 
     135      CASE ( np_FCT )                                    ! FCT      : 2nd / 4th order 
     136         CALL tra_adv_fct    ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra, nn_fct_h, nn_fct_v ) 
     137      CASE ( np_FCT_zts )                                ! 2nd order FCT with vertical time-splitting 
     138         CALL tra_adv_fct_zts( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra        , nn_fct_zts ) 
     139      CASE ( np_MUS )                                    ! MUSCL 
     140         CALL tra_adv_mus    ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra        , ln_mus_ups )  
     141      CASE ( np_UBS )                                    ! UBS 
     142         CALL tra_adv_ubs    ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra        , nn_ubs_v   ) 
     143      CASE ( np_QCK )                                    ! QUICKEST 
     144         CALL tra_adv_qck    ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra                     ) 
    123145      ! 
    124146      END SELECT 
     
    146168      INTEGER ::  ios                 ! Local integer output status for namelist read 
    147169      !! 
    148       NAMELIST/namtrc_adv/ ln_trcadv_cen2 , ln_trcadv_tvd   ,    & 
    149          &                 ln_trcadv_muscl, ln_trcadv_muscl2,    & 
    150          &                 ln_trcadv_ubs  , ln_trcadv_qck, ln_trcadv_msc_ups 
     170      NAMELIST/namtrc_adv/ ln_trcadv_cen, nn_cen_h, nn_cen_v,               &   ! CEN 
     171         &                 ln_trcadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts,   &   ! FCT 
     172         &                 ln_trcadv_mus,                     ln_mus_ups,   &   ! MUSCL 
     173         &                 ln_trcadv_ubs,           nn_ubs_v,               &   ! UBS 
     174         &                 ln_trcadv_qck                                        ! QCK 
    151175      !!---------------------------------------------------------------------- 
    152176      ! 
     
    165189         WRITE(numout,*) '~~~~~~~~~~~' 
    166190         WRITE(numout,*) '   Namelist namtrc_adv : chose a advection scheme for tracers' 
    167          WRITE(numout,*) '      2nd order advection scheme     ln_trcadv_cen2   = ', ln_trcadv_cen2 
    168          WRITE(numout,*) '      TVD advection scheme           ln_trcadv_tvd    = ', ln_trcadv_tvd 
    169          WRITE(numout,*) '      MUSCL  advection scheme        ln_trcadv_muscl  = ', ln_trcadv_muscl 
    170          WRITE(numout,*) '      MUSCL2 advection scheme        ln_trcadv_muscl2 = ', ln_trcadv_muscl2 
    171          WRITE(numout,*) '      UBS    advection scheme        ln_trcadv_ubs    = ', ln_trcadv_ubs 
    172          WRITE(numout,*) '      QUICKEST advection scheme      ln_trcadv_qck    = ', ln_trcadv_qck 
     191         WRITE(numout,*) '      centered scheme                           ln_trcadv_cen = ', ln_trcadv_cen 
     192         WRITE(numout,*) '            horizontal 2nd/4th order               nn_cen_h   = ', nn_fct_h 
     193         WRITE(numout,*) '            vertical   2nd/4th order               nn_cen_v   = ', nn_fct_v 
     194         WRITE(numout,*) '      Flux Corrected Transport scheme           ln_trcadv_fct = ', ln_trcadv_fct 
     195         WRITE(numout,*) '            horizontal 2nd/4th order               nn_fct_h   = ', nn_fct_h 
     196         WRITE(numout,*) '            vertical   2nd/4th order               nn_fct_v   = ', nn_fct_v 
     197         WRITE(numout,*) '            2nd order + vertical sub-timestepping  nn_fct_zts = ', nn_fct_zts 
     198         WRITE(numout,*) '      MUSCL scheme                              ln_trcadv_mus = ', ln_trcadv_mus 
     199         WRITE(numout,*) '            + upstream scheme near river mouths    ln_mus_ups = ', ln_mus_ups 
     200         WRITE(numout,*) '      UBS scheme                                ln_trcadv_ubs = ', ln_trcadv_ubs 
     201         WRITE(numout,*) '            vertical   2nd/4th order               nn_ubs_v   = ', nn_ubs_v 
     202         WRITE(numout,*) '      QUICKEST scheme                           ln_trcadv_qck = ', ln_trcadv_qck 
    173203      ENDIF 
    174204      ! 
    175205 
    176206      ioptio = 0                      ! Parameter control 
    177       IF( ln_trcadv_cen2   )   ioptio = ioptio + 1 
    178       IF( ln_trcadv_tvd    )   ioptio = ioptio + 1 
    179       IF( ln_trcadv_muscl  )   ioptio = ioptio + 1 
    180       IF( ln_trcadv_muscl2 )   ioptio = ioptio + 1 
    181       IF( ln_trcadv_ubs    )   ioptio = ioptio + 1 
    182       IF( ln_trcadv_qck    )   ioptio = ioptio + 1 
    183       ! 
     207      IF( ln_trcadv_cen                      )   nadv = np_CEN 
     208      IF( ln_trcadv_fct                      )   nadv = np_FCT 
     209      IF( ln_trcadv_fct .AND. nn_fct_zts > 0 )   nadv = np_FCT_zts 
     210      IF( ln_trcadv_mus                      )   nadv = np_MUS 
     211      IF( ln_trcadv_ubs                      )   nadv = np_UBS 
     212      IF( ln_trcadv_qck                      )   nadv = np_QCK 
     213      ! 
     214      IF( ioptio == 0 ) THEN 
     215         nadv = np_NO_adv 
     216         CALL ctl_warn( 'trc_adv_init: You are running without tracer advection.' ) 
     217      ENDIF 
    184218      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtrc_adv' ) 
    185219      ! 
     220      IF( ln_trcadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 )   & 
     221                        .AND. ( nn_cen_v /= 2 .AND. nn_cen_v /= 4 )   ) THEN 
     222        CALL ctl_stop( 'trc_adv_init: CEN scheme, choose 2nd or 4th order' ) 
     223      ENDIF 
     224      IF( ln_trcadv_fct .AND. ( nn_fct_h /= 2 .AND. nn_fct_h /= 4 )   & 
     225                        .AND. ( nn_fct_v /= 2 .AND. nn_fct_v /= 4 )   ) THEN 
     226        CALL ctl_stop( 'trc_adv_init: FCT scheme, choose 2nd or 4th order' ) 
     227      ENDIF 
     228      IF( ln_trcadv_fct .AND. nn_fct_zts > 0 ) THEN 
     229         IF( nn_fct_h == 4 ) THEN 
     230            nn_fct_h = 2 
     231            CALL ctl_stop( 'trc_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' ) 
     232         ENDIF 
     233         IF( lk_vvl ) THEN 
     234            CALL ctl_stop( 'trc_adv_init: vertical sub-timestepping not allow in non-linear free surface' ) 
     235         ENDIF 
     236         IF( nn_fct_zts == 1 )   CALL ctl_warn( 'trc_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' ) 
     237      ENDIF 
     238      IF( ln_trcadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 )   ) THEN 
     239        CALL ctl_stop( 'trc_adv_init: UBS scheme, choose 2nd or 4th order' ) 
     240      ENDIF 
     241      IF( ln_trcadv_ubs .AND. nn_ubs_v == 4 ) THEN 
     242         CALL ctl_warn( 'trc_adv_init: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' ) 
     243      ENDIF 
     244      IF( ln_isfcav ) THEN                                                       ! ice-shelf cavities 
     245         IF(  ln_trcadv_cen .AND. nn_cen_v /= 4    .OR.   &                            ! NO 4th order with ISF 
     246            & ln_trcadv_fct .AND. nn_fct_v /= 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 
     247      ENDIF 
     248      ! 
    186249      !                              ! Set nadv 
    187       IF( ln_trcadv_cen2   )   nadv =  1 
    188       IF( ln_trcadv_tvd    )   nadv =  2 
    189       IF( ln_trcadv_muscl  )   nadv =  3 
    190       IF( ln_trcadv_muscl2 )   nadv =  4 
    191       IF( ln_trcadv_ubs    )   nadv =  5 
    192       IF( ln_trcadv_qck    )   nadv =  6 
     250      IF( ln_trcadv_cen  )   nadv = np_CEN 
     251      IF( ln_trcadv_fct  )   nadv = np_FCT 
     252      IF( nn_fct_zts > 0 )   nadv = np_FCT_zts 
     253      IF( ln_trcadv_mus  )   nadv = np_MUS 
     254      IF( ln_trcadv_ubs  )   nadv = np_UBS 
     255      IF( ln_trcadv_qck  )   nadv = np_QCK 
    193256      ! 
    194257      IF(lwp) THEN                   ! Print the choice 
    195258         WRITE(numout,*) 
    196          IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used' 
    197          IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used' 
    198          IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used' 
    199          IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    scheme is used' 
    200          IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used' 
    201          IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used' 
     259         IF( nadv == np_NO_adv  )   WRITE(numout,*) '         NO passive tracer advection' 
     260         IF( nadv == np_CEN     )   WRITE(numout,*) '         CEN      scheme is used. Horizontal order: ', nn_cen_h,   & 
     261            &                                                                        ' Vertical   order: ', nn_cen_v 
     262         IF( nadv == np_FCT     )   WRITE(numout,*) '         FCT      scheme is used. Horizontal order: ', nn_fct_h,   & 
     263            &                                                                       ' Vertical   order: ', nn_fct_v 
     264         IF( nadv == np_FCT_zts )   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 
     265         IF( nadv == np_MUS     )   WRITE(numout,*) '         MUSCL    scheme is used' 
     266         IF( nadv == np_UBS     )   WRITE(numout,*) '         UBS      scheme is used' 
     267         IF( nadv == np_QCK     )   WRITE(numout,*) '         QUICKEST scheme is used' 
    202268      ENDIF 
    203269      ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90

    r5766 r5770  
    7373      IF( nn_timing == 1 )   CALL timing_start('trc_ldf') 
    7474      ! 
    75        
    7675      IF( l_trdtrc )  THEN 
    7776         CALL wrk_alloc( jpi,jpj,jpk,jptra,   ztrtrd ) 
    7877         ztrtrd(:,:,:,:)  = tra(:,:,:,:) 
    7978      ENDIF 
    80  
    81       !                                        ! set the lateral diffusivity coef. for passive tracer       
     79      ! 
     80      !                                        !* set the lateral diffusivity coef. for passive tracer       
    8281      CALL wrk_alloc( jpi,jpj,jpk,   zahu, zahv ) 
    8382      zahu(:,:,:) = rldf * ahtu(:,:,:) 
     
    8685      SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend 
    8786      ! 
    88       CASE ( n_lap   )                                ! iso-level laplacian 
     87      CASE ( np_lap   )                               ! iso-level laplacian 
    8988         CALL tra_ldf_lap  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb,      tra, jptra,  1   ) 
    9089         ! 
    91       CASE ( n_lap_i )                                ! laplacian : standard iso-neutral operator (Madec) 
     90      CASE ( np_lap_i )                               ! laplacian : standard iso-neutral operator (Madec) 
    9291         CALL tra_ldf_iso  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,  1   ) 
    9392         ! 
    94       CASE ( n_lap_it )                               ! laplacian : triad iso-neutral operator (griffies) 
     93      CASE ( np_lap_it )                              ! laplacian : triad iso-neutral operator (griffies) 
    9594         CALL tra_ldf_triad( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,  1   ) 
    9695         ! 
    97       CASE ( n_blp , n_blp_i , n_blp_it )             ! bilaplacian: all operator (iso-level, -neutral) 
     96      CASE ( np_blp , np_blp_i , np_blp_it )          ! bilaplacian: all operator (iso-level, -neutral) 
    9897         CALL tra_ldf_blp  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb     , tra, jptra, nldf ) 
    9998         ! 
     
    119118   END SUBROUTINE trc_ldf 
    120119 
    121 !!gm ldf_ctl should be called in trcini  so that l_ldfslp=T  cause the slope init and calculation 
    122120 
    123121   SUBROUTINE trc_ldf_ini 
     
    128126      !! 
    129127      !! ** Method  :   set nldf from the namtra_ldf logicals 
    130       !!      nldf == -2   No lateral diffusion 
    131128      !!      nldf ==  0   laplacian operator 
    132129      !!      nldf ==  1   Rotated laplacian operator 
     
    134131      !!      nldf ==  3   Rotated bilaplacian 
    135132      !!---------------------------------------------------------------------- 
    136       INTEGER ::   ioptio, ierr         ! temporary integers 
    137       INTEGER ::  ios                 ! Local integer output status for namelist read 
    138       ! 
    139       NAMELIST/namtrc_ldf/ ln_trcldf_lap, ln_trcldf_blp, & 
     133      INTEGER ::   ioptio, ierr   ! temporary integers 
     134      INTEGER ::   ios            ! Local integer output status for namelist read 
     135      ! 
     136      NAMELIST/namtrc_ldf/ ln_trcldf_lap, ln_trcldf_blp,                                  & 
    140137         &                 ln_trcldf_lev, ln_trcldf_hor, ln_trcldf_iso, ln_trcldf_triad,  & 
    141138         &                 rn_ahtrc_0   , rn_bhtrc_0 
     
    172169      IF( ln_trcldf_lap )   ioptio = ioptio + 1 
    173170      IF( ln_trcldf_blp )   ioptio = ioptio + 1 
    174       IF( ioptio >  1 )   CALL ctl_stop( 'trc_ldf_ctl: use ONE or NONE of the 2 lap/bilap operator type on tracer' ) 
    175       IF( ioptio == 0 )   nldf = n_no_ldf   ! No lateral diffusion 
     171      IF( ioptio >  1   )   CALL ctl_stop( 'trc_ldf_ctl: use ONE or NONE of the 2 lap/bilap operator type on tracer' ) 
     172      IF( ioptio == 0   )   nldf = np_no_ldf   ! No lateral diffusion 
    176173       
    177174      IF( ln_trcldf_lap .AND. ln_trcldf_blp )   CALL ctl_stop( 'trc_ldf_ctl: bilaplacian should be used on both TRC and TRA' ) 
     
    189186      IF( ln_trcldf_lap ) THEN      !==  laplacian operator  ==! 
    190187         IF ( ln_zco ) THEN                ! z-coordinate 
    191             IF ( ln_trcldf_lev   )   nldf = n_lap     ! iso-level = horizontal (no rotation) 
    192             IF ( ln_trcldf_hor   )   nldf = n_lap     ! iso-level = horizontal (no rotation) 
    193             IF ( ln_trcldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard  (   rotation) 
    194             IF ( ln_trcldf_triad )   nldf = n_lap_it  ! iso-neutral: triad     (   rotation) 
     188            IF ( ln_trcldf_lev   )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     189            IF ( ln_trcldf_hor   )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     190            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation) 
     191            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation) 
    195192         ENDIF 
    196193         IF ( ln_zps ) THEN             ! z-coordinate with partial step 
    197194            IF ( ln_trcldf_lev   )   ierr = 1         ! iso-level not allowed  
    198             IF ( ln_trcldf_hor   )   nldf = n_lap     ! horizontal (no rotation) 
    199             IF ( ln_trcldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation) 
    200             IF ( ln_trcldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation) 
     195            IF ( ln_trcldf_hor   )   nldf = np_lap     ! horizontal (no rotation) 
     196            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard (rotation) 
     197            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad    (rotation) 
    201198         ENDIF 
    202199         IF ( ln_sco ) THEN             ! s-coordinate 
    203             IF ( ln_trcldf_lev   )   nldf = n_lap     ! iso-level  (no rotation) 
    204             IF ( ln_trcldf_hor   )   nldf = n_lap_it  ! horizontal (   rotation)       !!gm   a checker.... 
    205             IF ( ln_trcldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation) 
    206             IF ( ln_trcldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation) 
     200            IF ( ln_trcldf_lev   )   nldf = np_lap     ! iso-level  (no rotation) 
     201            IF ( ln_trcldf_hor   )   nldf = np_lap_it  ! horizontal (   rotation)       !!gm   a checker.... 
     202            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard (rotation) 
     203            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad    (rotation) 
    207204         ENDIF 
    208205         !                                ! diffusivity ratio: passive / active tracers  
     
    217214         ENDIF 
    218215      ENDIF 
    219  
     216      ! 
    220217      IF( ln_trcldf_blp ) THEN      !==  bilaplacian operator  ==! 
    221218         IF ( ln_zco ) THEN                ! z-coordinate 
    222             IF ( ln_trcldf_lev   )   nldf = n_blp     ! iso-level = horizontal (no rotation) 
    223             IF ( ln_trcldf_hor   )   nldf = n_blp     ! iso-level = horizontal (no rotation) 
    224             IF ( ln_trcldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
    225             IF ( ln_trcldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     219            IF ( ln_trcldf_lev   )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
     220            IF ( ln_trcldf_hor   )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
     221            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation) 
     222            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation) 
    226223         ENDIF 
    227224         IF ( ln_zps ) THEN             ! z-coordinate with partial step 
    228225            IF ( ln_trcldf_lev   )   ierr = 1         ! iso-level not allowed  
    229             IF ( ln_trcldf_hor   )   nldf = n_blp     ! horizontal (no rotation) 
    230             IF ( ln_trcldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
    231             IF ( ln_trcldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     226            IF ( ln_trcldf_hor   )   nldf = np_blp     ! horizontal (no rotation) 
     227            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation) 
     228            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation) 
    232229         ENDIF 
    233230         IF ( ln_sco ) THEN             ! s-coordinate 
    234             IF ( ln_trcldf_lev   )   nldf = n_blp     ! iso-level  (no rotation) 
    235             IF ( ln_trcldf_hor   )   nldf = n_blp_it  ! horizontal (   rotation)       !!gm   a checker.... 
    236             IF ( ln_trcldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
    237             IF ( ln_trcldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     231            IF ( ln_trcldf_lev   )   nldf = np_blp     ! iso-level  (no rotation) 
     232            IF ( ln_trcldf_hor   )   nldf = np_blp_it  ! horizontal (   rotation)       !!gm   a checker.... 
     233            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation) 
     234            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation) 
    238235         ENDIF 
    239236         !                                ! diffusivity ratio: passive / active tracers  
     
    248245         ENDIF 
    249246      ENDIF 
    250  
     247      ! 
    251248      IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' ) 
    252249      IF( ln_ldfeiv .AND. .NOT.ln_trcldf_iso )   & 
     
    256253         IF( .NOT.l_ldfslp )   CALL ctl_stop( '          the rotation of the diffusive tensor require l_ldfslp' ) 
    257254      ENDIF 
    258  
     255      ! 
    259256      IF(lwp) THEN 
    260257         WRITE(numout,*) 
    261          IF( nldf == n_no_ldf )   WRITE(numout,*) '          NO lateral diffusion' 
    262          IF( nldf == n_lap    )   WRITE(numout,*) '          laplacian iso-level operator' 
    263          IF( nldf == n_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)' 
    264          IF( nldf == n_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)' 
    265          IF( nldf == n_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator' 
    266          IF( nldf == n_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)' 
    267          IF( nldf == n_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)' 
     258         IF( nldf == np_no_ldf )   WRITE(numout,*) '          NO lateral diffusion' 
     259         IF( nldf == np_lap    )   WRITE(numout,*) '          laplacian iso-level operator' 
     260         IF( nldf == np_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)' 
     261         IF( nldf == np_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)' 
     262         IF( nldf == np_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator' 
     263         IF( nldf == np_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)' 
     264         IF( nldf == np_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)' 
    268265      ENDIF 
    269266      ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90

    r5766 r5770  
    1111   !!   'key_top'                                                TOP models 
    1212   !!---------------------------------------------------------------------- 
    13    !!   trc_zdf     : update the tracer trend with the lateral diffusion 
    14    !!   trc_zdf_ini : initialization, namelist read, and parameters control 
     13   !!   trc_zdf      : update the tracer trend with the lateral diffusion 
     14   !!   trc_zdf_ini  : initialization, namelist read, and parameters control 
    1515   !!---------------------------------------------------------------------- 
    16    USE oce_trc         ! ocean dynamics and active tracers 
    17    USE trc             ! ocean passive tracers variables 
    18    USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp     routine) 
    19    USE trazdf_imp      ! vertical diffusion: implicit (tra_zdf_imp     routine) 
    20    USE trcldf    
    21    USE trd_oce 
    22    USE trdtra 
    23    USE prtctl_trc      ! Print control 
     16   USE trc           ! ocean passive tracers variables 
     17   USE oce_trc       ! ocean dynamics and active tracers 
     18   USE trd_oce       ! trends: ocean variables 
     19   USE trazdf_exp    ! vertical diffusion: explicit (tra_zdf_exp     routine) 
     20   USE trazdf_imp    ! vertical diffusion: implicit (tra_zdf_imp     routine) 
     21   USE trcldf        ! passive tracers: lateral diffusion 
     22   USE trdtra        ! trends manager: tracers  
     23   USE prtctl_trc    ! Print control 
    2424 
    2525   IMPLICIT NONE 
    2626   PRIVATE 
    2727 
    28    PUBLIC   trc_zdf          ! called by step.F90  
    29    PUBLIC   trc_zdf_alloc    ! called by nemogcm.F90  
    30    PUBLIC   trc_zdf_ini    ! called by nemogcm.F90  
    31    !                                        !!: ** Vertical diffusion 
    32    !                                        (nam_trczdf) ** 
     28   PUBLIC   trc_zdf         ! called by step.F90  
     29   PUBLIC   trc_zdf_ini     ! called by nemogcm.F90  
     30   PUBLIC   trc_zdf_alloc   ! called by nemogcm.F90  
     31    
     32   !                                        !!** Vertical diffusion (nam_trczdf) ** 
    3333   LOGICAL , PUBLIC ::   ln_trczdf_exp       !: explicit vertical diffusion scheme flag 
    3434   INTEGER , PUBLIC ::   nn_trczdf_exp       !: number of sub-time step (explicit time stepping) 
     
    4444#  include "vectopt_loop_substitute.h90" 
    4545   !!---------------------------------------------------------------------- 
    46    !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     46   !! NEMO/TOP 3.7 , NEMO Consortium (2015) 
    4747   !! $Id$  
    4848   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    5959      ! 
    6060   END FUNCTION trc_zdf_alloc 
     61 
    6162 
    6263   SUBROUTINE trc_zdf( kt ) 
     
    8788 
    8889      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend 
    89       CASE ( -1 )                                       ! esopa: test all possibility with control print 
    90          CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )  
    91          WRITE(charout, FMT="('zdf1 ')") ;  CALL prt_ctl_trc_info(charout) 
    92                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 
    93          CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dt,                trb, tra, jptra )  
    94          WRITE(charout, FMT="('zdf2 ')") ;  CALL prt_ctl_trc_info(charout) 
    95                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 
    9690      CASE ( 0 ) ;  CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme  
    9791      CASE ( 1 ) ;  CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dt,                trb, tra, jptra )    !   implicit scheme           
    98  
    9992      END SELECT 
    10093 
     
    118111   END SUBROUTINE trc_zdf 
    119112 
     113 
    120114   SUBROUTINE trc_zdf_ini 
    121115      !!---------------------------------------------------------------------- 
     
    132126      !!---------------------------------------------------------------------- 
    133127      INTEGER ::  ios                 ! Local integer output status for namelist read 
    134       ! 
     128      !! 
    135129      NAMELIST/namtrc_zdf/ ln_trczdf_exp  , nn_trczdf_exp 
    136130      !!---------------------------------------------------------------------- 
    137  
    138       !                                ! Vertical mixing 
     131      ! 
    139132      REWIND( numnat_ref )             ! namtrc_zdf in reference namelist  
    140133      READ  ( numnat_ref, namtrc_zdf, IOSTAT = ios, ERR = 905) 
    141134905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_zdf in reference namelist', lwp ) 
    142  
     135      ! 
    143136      REWIND( numnat_cfg )             ! namtrc_zdf in configuration namelist  
    144137      READ  ( numnat_cfg, namtrc_zdf, IOSTAT = ios, ERR = 906 ) 
    145138906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_zdf in configuration namelist', lwp ) 
    146139      IF(lwm) WRITE ( numont, namtrc_zdf ) 
    147  
    148       IF(lwp) THEN                     !   ! Control print 
     140      ! 
     141      IF(lwp) THEN                     ! Control print 
    149142         WRITE(numout,*) 
    150143         WRITE(numout,*) '   Namelist namtrc_zdf : set vertical diffusion  parameters' 
     
    153146      ENDIF 
    154147 
    155       !  Define the vertical tracer physics scheme 
    156       IF( ln_trczdf_exp ) THEN  ;  nzdf = 0     ! explicit scheme 
    157       ELSE                      ;  nzdf = 1     ! implicit scheme 
     148      !                                ! Define the vertical tracer physics scheme 
     149      IF( ln_trczdf_exp ) THEN   ;   nzdf = 0     ! explicit scheme 
     150      ELSE                       ;   nzdf = 1     ! implicit scheme 
    158151      ENDIF 
    159152 
    160       ! Force implicit schemes 
    161       IF( ln_trcldf_iso                               )   nzdf = 1      ! iso-neutral lateral physics 
    162       IF( ln_trcldf_hor .AND. ln_sco                  )   nzdf = 1      ! horizontal lateral physics in s-coordinate 
     153      !                                ! Force implicit schemes 
     154      IF( ln_trcldf_iso              )   nzdf = 1      ! iso-neutral lateral physics 
     155      IF( ln_trcldf_hor .AND. ln_sco )   nzdf = 1      ! horizontal lateral physics in s-coordinate 
    163156#if defined key_zdftke || defined key_zdfgls  
    164                                                           nzdf = 1      ! TKE or GLS physics        
     157                                         nzdf = 1      ! TKE or GLS physics        
    165158#endif 
    166159      IF( ln_trczdf_exp .AND. nzdf == 1 )  &  
     
    175168         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme' 
    176169      ENDIF 
    177  
     170      ! 
    178171   END SUBROUTINE trc_zdf_ini 
     172    
    179173#else 
    180174   !!---------------------------------------------------------------------- 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r5766 r5770  
    5353      !!                or read data or analytical formulation 
    5454      !!--------------------------------------------------------------------- 
    55       !!--------------------------------------------------------------------- 
    5655      ! 
    5756      IF( nn_timing == 1 )   CALL timing_start('trc_init') 
     
    9291      ! 
    9392   END SUBROUTINE trc_init 
     93 
    9494 
    9595   SUBROUTINE trc_ini_ctl 
     
    103103      l_trcdm2dc = ln_dm2dc .OR. ( ln_cpl .AND. ncpl_qsr_freq /= 1 ) 
    104104      l_trcdm2dc = l_trcdm2dc  .AND. .NOT. lk_offline 
    105       IF( l_trcdm2dc .AND. lwp ) & 
    106          &   CALL ctl_warn(' Coupling with passive tracers and used of diurnal cycle. & 
    107          & Computation of a daily mean shortwave for some biogeochemical models) ') 
    108  
    109       ! 
    110       IF( nn_cla == 1 )   & 
    111          &  CALL ctl_stop( ' Cross Land Advection not yet implemented with passive tracer ; nn_cla must be 0' ) 
     105      IF( l_trcdm2dc .AND. lwp )   CALL ctl_warn( 'Coupling with passive tracers and used of diurnal cycle.',   & 
     106         &                           'Computation of a daily mean shortwave for some biogeochemical models ' ) 
    112107      ! 
    113108   END SUBROUTINE trc_ini_ctl 
     109 
    114110 
    115111   SUBROUTINE trc_ini_inv 
     
    157153   END SUBROUTINE trc_ini_inv 
    158154 
     155 
    159156   SUBROUTINE trc_ini_sms 
    160157      !!---------------------------------------------------------------------- 
     
    196193      ! 
    197194   END SUBROUTINE trc_ini_trp 
     195 
    198196 
    199197   SUBROUTINE trc_ini_state 
     
    248246   END SUBROUTINE trc_ini_state 
    249247 
     248 
    250249   SUBROUTINE top_alloc 
    251250      !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.