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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6772 for branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Ignore:
Timestamp:
2016-07-01T18:02:45+02:00 (8 years ago)
Author:
cbricaud
Message:

clean in coarsening branch

Location:
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2_crs.F90

    r5601 r6772  
    6161   ! 
    6262   PUBLIC   eos_crs        ! called by step, istate, tranpc and zpsgrd modules 
    63    PUBLIC   bn2_crs        ! called by step module 
    6463   PUBLIC   eos_rab_crs    ! called by ldfslp, zdfddm, trabbl 
    6564   PUBLIC   eos_init_crs   ! called by istate module 
    … …  
    392391               DO ji = 1, jpi_crs 
    393392                  ! 
    394                   zh  = gdept_crs(ji,jj,jk) * r1_Z0                                ! depth 
     393                  zh  = fsdept_crs(ji,jj,jk) * r1_Z0                                ! depth 
    395394                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    396395                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    … …  
    450449                  zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    451450                  zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    452                   zh  = gdept_crs(ji,jj,jk)                 ! depth in meters at t-point 
     451                  zh  = fsdept_crs(ji,jj,jk)                 ! depth in meters at t-point 
    453452                  ztm = tmask_crs(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
    454453                  ! 
    … …  
    689688      ! 
    690689   END SUBROUTINE rab_crs_0d 
    691  
    692  
    693    SUBROUTINE bn2_crs( pts, pab, pn2 ) 
    694       !!---------------------------------------------------------------------- 
    695       !!                  ***  ROUTINE bn2  *** 
    696       !! 
    697       !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the  
    698       !!                time-step of the input arguments 
    699       !! 
    700       !! ** Method  :   pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w 
    701       !!      where alpha and beta are given in pab, and computed on T-points. 
    702       !!      N.B. N^2 is set one for all to zero at jk=1 in istate module. 
    703       !! 
    704       !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point  
    705       !! 
    706       !!---------------------------------------------------------------------- 
    707       REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu] 
    708       REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celcius-1,psu-1] 
    709       REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2] 
    710       ! 
    711       INTEGER  ::   ji, jj, jk      ! dummy loop indices 
    712       REAL(wp) ::   zaw, zbw, zrw   ! local scalars 
    713       !!---------------------------------------------------------------------- 
    714       ! 
    715       pn2(:,:,:)=0._wp 
    716  
    717       IF( nn_timing == 1 ) CALL timing_start('bn2') 
    718       ! 
    719       DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 ) 
    720          DO jj = 1, jpj_crs          ! surface and bottom value set to zero one for all in istate.F90 
    721             DO ji = 1, jpi_crs 
    722                !zrw =   ( gdepw_crs(ji,jj,jk  ) - gdept_crs(ji,jj,jk) )   & 
    723                !   &  / ( gdept_crs(ji,jj,jk-1) - gdept_crs(ji,jj,jk) )  
    724                zrw =   gdepw_crs(ji,jj,jk  ) - gdept_crs(ji,jj,jk)     
    725                !?IF( gdept_crs(ji,jj,jk-1) - gdept_crs(ji,jj,jk) .NE. 0._wp )THEN 
    726                IF( gdept_crs(ji,jj,jk-1) - gdept_crs(ji,jj,jk) .LT. 0._wp )THEN 
    727                   zrw = zrw  / ( gdept_crs(ji,jj,jk-1) - gdept_crs(ji,jj,jk) )  
    728                ELSE 
    729                   zrw = 0._wp 
    730                ENDIF 
    731                ! 
    732                zaw = pab(ji,jj,jk,jp_tem) * (1._wp - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
    733                zbw = pab(ji,jj,jk,jp_sal) * (1._wp - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
    734                ! 
    735                IF( e3w_max_crs(ji,jj,jk) .NE. 0._wp ) THEN 
    736                   pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
    737                      &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
    738                   &                    * tmask_crs(ji,jj,jk)  / e3w_max_crs(ji,jj,jk) 
    739                ENDIF 
    740             END DO 
    741          END DO 
    742       END DO 
    743       ! 
    744       IF( nn_timing == 1 )   CALL timing_stop('bn2') 
    745       ! 
    746    END SUBROUTINE bn2_crs 
    747690 
    748691   SUBROUTINE eos_init_crs 

  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd_crs.F90

    r6101 r6772  
    9191      !!---------------------------------------------------------------------- 
    9292      ! 
     93 
    9394      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd') 
    9495      ! 
    … …  
    126127         ! upstream tracer flux in the i and j direction 
    127128         DO jk = 1, jpkm1 
    128             DO jj = 1, jpjm1 
    129                DO ji = 1, fs_jpim1   ! vector opt. 
     129            DO jj = 2, jpj_crs-1 
     130               DO ji = 2, jpi_crs-1 
    130131                  ! upstream scheme 
    131132                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
    … …  
    138139            END DO 
    139140         END DO 
     141         CALL crs_lbc_lnk( zwx, 'U', -1._wp )   
     142         CALL crs_lbc_lnk( zwy, 'V', -1._wp )   
    140143         ! upstream tracer flux in the k direction 
    141144         ! Surface value 
    142145         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable 
    143          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface  
     146         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) !cbr * ptb(:,:,1,jn)   ! linear free surface  
    144147         ENDIF 
    145148         ! Interior value 
    146149         DO jk = 2, jpkm1 
    147             DO jj = 1, jpj 
    148                DO ji = 1, jpi 
     150            DO jj = 2,  jpj_crs-1 
     151               DO ji = nldi_crs, nlei_crs 
    149152                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    150153                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    … …  
    153156            END DO 
    154157         END DO 
     158         CALL crs_lbc_lnk( zwz, 'T', 1. )   
     159 
    155160         ! total advective trend 
    156161         DO jk = 1, jpkm1 
    157162            z2dtt = p2dt(jk) 
    158             DO jj = 2, jpjm1 
    159                DO ji = fs_2, fs_jpim1   ! vector opt. 
     163            DO jj = 2, jpj_crs-1 
     164               DO ji = 2, jpi_crs-1 
    160165                  zbtr = r1_bt_crs(ji,jj,jk)  
    161166                  ! total intermediate advective trends 
    … …  
    163168                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    164169                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
    165                   ! update and guess with monotonic sheme 
     170 
    166171                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra 
    167172                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask_crs(ji,jj,jk) 
    … …  
    169174            END DO 
    170175         END DO 
     176  
    171177         !                             ! Lateral boundary conditions on zwi  (unchanged sign) 
    172178         CALL crs_lbc_lnk( zwi, 'T', 1. )   
    … …  
    187193         ! antidiffusive flux on i and j 
    188194         DO jk = 1, jpkm1 
    189             DO jj = 1, jpjm1 
    190                DO ji = 1, fs_jpim1   ! vector opt. 
     195            DO jj = 2, jpj_crs-1 
     196               DO ji = 2, jpi_crs-1 
    191197                  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) 
    192198                  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) 
    … …  
    198204         ! 
    199205         DO jk = 2, jpkm1          ! Interior value 
    200             DO jj = 1, jpj 
    201                DO ji = 1, jpi 
     206            DO jj = 2, jpj_crs-1 
     207               DO ji = 2, jpi_crs-1 
    202208                  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) 
    203209               END DO 
    204210            END DO 
    205          END DO 
     211        END DO 
    206212         CALL crs_lbc_lnk( zwx, 'U', -1. )   ;   CALL crs_lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    207213         CALL crs_lbc_lnk( zwz, 'W',  1. ) 
    … …  
    214220         ! ------------------------------------ 
    215221         DO jk = 1, jpkm1 
    216             DO jj = 2, jpjm1 
    217                DO ji = fs_2, fs_jpim1   ! vector opt.   
     222            DO jj = 2, jpj_crs-1 
     223               DO ji = 2, jpi_crs-1 
    218224                  zbtr = r1_bt_crs(ji,jj,jk) 
    219225                  ! total advective trends 
    … …  
    247253      END DO 
    248254      ! 
     255 
    249256                   CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz , zwx, zwy ) 
    250257      IF( l_trd )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
    … …  
    302309         ikm1 = MAX(jk-1,1) 
    303310         z2dtt = p2dt(jk) 
    304          DO jj = 2, jpjm1 
    305             DO ji = fs_2, fs_jpim1   ! vector opt. 
     311         DO jj = 2, jpj_crs-1 
     312            DO ji = 2, jpi_crs-1 
    306313 
    307314               ! search maximum in neighbourhood 
    … …  
    339346      ! ---------------------------------------- 
    340347      DO jk = 1, jpkm1 
    341          DO jj = 2, jpjm1 
    342             DO ji = fs_2, fs_jpim1   ! vector opt. 
     348         DO jj = 2, jpj_crs-1 
     349            DO ji = 2, jpi_crs-1 
    343350               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
    344351               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 

  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_crs.F90

    r6101 r6772  
    1010   !!            3.3  !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
    1111   !!---------------------------------------------------------------------- 
    12 #if   defined key_ldfslp   ||   defined key_esopa 
     12#if   ( defined key_ldfslp   ||   defined key_esopa ) && defined key_crs 
    1313   !!---------------------------------------------------------------------- 
    1414   !!   'key_ldfslp'               slope of the lateral diffusive direction 
    … …  
    1919   !!                  the isopycnal or geopotential s-coord. operator  
    2020   !!---------------------------------------------------------------------- 
    21 !   USE oce             ! ocean dynamics and active tracers 
    22 !   USE dom_oce         ! ocean space and time domain 
    23 !   USE trc_oce         ! share passive tracers/Ocean variables 
    24 !   USE zdf_oce         ! ocean vertical physics 
    25 !   USE ldftra_oce      ! ocean active tracers: lateral physics 
    26 !   USE ldfslp          ! iso-neutral slopes 
    2721   USE ldfslp_crs          ! iso-neutral slopes 
    2822   USE diaptr          ! poleward transport diagnostics 
    … …  
    3529   USE wrk_nemo        ! Memory Allocation 
    3630   USE timing          ! Timing 
    37 !   USE crs 
    3831   USE oce_trc 
    3932   USE iom, ONLY : iom_put,iom_swap 
    … …  
    113106      REAL(wp)                         ::   zztmp               ! local scalar 
    114107#endif 
     108      REAL(wp)                         ::   zmin,zmax 
    115109      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zdkt, zdk1t, z2d 
    116110      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdit, zdjt, ztfw , zftu,  zftv  
    … …  
    188182                     &             + tmask_crs(ji,jj+1,jk+1) + tmask_crs(ji,jj,jk  ), 1. ) 
    189183                  ! 
    190                   zcof1 = - fsahtu(ji,jj,jk) * e2e3u_msk(ji,jj,jk) * uslp_crs(ji,jj,jk) * zmsku / MAX( 1._wp , e3u_max_crs(ji,jj,jk)) 
    191                   zcof2 = - fsahtv(ji,jj,jk) * e1e3v_msk(ji,jj,jk) * vslp_crs(ji,jj,jk) * zmskv / MAX( 1._wp , e3v_max_crs(ji,jj,jk)) 
     184                  zcof1 = - fsahtu(ji,jj,jk) * e2e3u_msk(ji,jj,jk) * uslp_crs(ji,jj,jk) * zmsku / MAX( 1._wp , fse3u_max_crs(ji,jj,jk)) 
     185                  zcof2 = - fsahtv(ji,jj,jk) * e1e3v_msk(ji,jj,jk) * vslp_crs(ji,jj,jk) * zmskv / MAX( 1._wp , fse3v_max_crs(ji,jj,jk)) 
    192186                  ! 
    193187                  zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk)   & 
    … …  
    199193               END DO 
    200194            END DO 
    201             CALL iom_swap( "nemo_crs"  ) 
    202             CALL iom_put( "zftu" , zftu ) 
    203             CALL iom_put( "zftv" , zftv ) 
    204             CALL iom_swap( "nemo" ) 
    205195 
    206196            ! II.4 Second derivative (divergence) and add to the general trend 

  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_crs.F90

    r6101 r6772  
    124124                     ikv = mbkv_crs(ji,jj) 
    125125                     IF( iku == jk ) THEN 
    126                         zabe1 = fsahtu(ji,jj,iku) * umask_crs(ji,jj,iku) * e1ur(ji,jj) * e3u_crs(ji,jj,iku) 
     126                        zabe1 = fsahtu(ji,jj,iku) * umask_crs(ji,jj,iku) * e1ur(ji,jj) * fse3u_crs(ji,jj,iku) 
    127127                        ztu(ji,jj,jk) = zabe1 * pgu(ji,jj,jn) 
    128128                     ENDIF 
    129129                     IF( ikv == jk ) THEN 
    130                         zabe2 = fsahtv(ji,jj,ikv) * vmask_crs(ji,jj,ikv) * e2vr(ji,jj) * e3v_crs(ji,jj,ikv) 
     130                        zabe2 = fsahtv(ji,jj,ikv) * vmask_crs(ji,jj,ikv) * e2vr(ji,jj) * fse3v_crs(ji,jj,ikv) 
    131131                        ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 
    132132                     ENDIF 

  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r5602 r6772  
    4949   USE agrif_opa_interp 
    5050#endif 
     51   USE crs 
    5152 
    5253   IMPLICIT NONE 
    … …  
    5657   PUBLIC   tra_nxt_fix   ! to be used in trcnxt 
    5758   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt 
     59   PUBLIC   tra_nxt_vvl_crs ! to be used in trcnxt 
    5860 
    5961   REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg 
    … …  
    349351   END SUBROUTINE tra_nxt_vvl 
    350352 
     353  SUBROUTINE tra_nxt_vvl_crs( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt ) 
     354      !!---------------------------------------------------------------------- 
     355      !!                   ***  ROUTINE tra_nxt_vvl  *** 
     356      !! 
     357      !! ** Purpose :   Time varying volume: apply the Asselin time filter   
     358      !!                and swap the tracer fields. 
     359      !!  
     360      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
     361      !!              - save in (ta,sa) a thickness weighted average over the three  
     362      !!             time levels which will be used to compute rdn and thus the semi- 
     363      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 
     364      !!              - swap tracer fields to prepare the next time_step. 
     365      !!                This can be summurized for tempearture as: 
     366      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T 
     367      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )    
     368      !!             ztm = 0                                                       otherwise 
     369      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
     370      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] ) 
     371      !!             tn  = ta  
     372      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call) 
     373      !! 
     374      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
     375      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     376      !!---------------------------------------------------------------------- 
     377      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index 
     378      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index 
     379      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step 
     380      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator) 
     381      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers 
     382      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields 
     383      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields 
     384      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend 
     385      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content 
     386      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content 
     387 
     388      !!      
     389      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical 
     390      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
     391      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar 
     392      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      - 
     393      !!---------------------------------------------------------------------- 
     394      !!---------------------------------------------------------------------- 
     395      ! 
     396      IF( kt == kit000 )  THEN 
     397         IF(lwp) WRITE(numout,*) 
     398         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype 
     399         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     400      ENDIF 
     401      ! 
     402      IF( cdtype == 'TRA' )  THEN 
     403         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg 
     404         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration 
     405         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs 
     406      ELSE 
     407         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg 
     408         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration 
     409         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs 
     410      ENDIF 
     411      ! 
     412      DO jn = 1, kjpt 
     413         DO jk = 1, jpkm1 
     414            zfact1 = atfp * p2dt(jk) 
     415            zfact2 = zfact1 / rau0 
     416            DO jj = 1, jpj 
     417               DO ji = 1, jpi 
     418                  ze3t_b = fse3t_b_crs(ji,jj,jk) 
     419                  ze3t_n = fse3t_n_crs(ji,jj,jk) 
     420                  ze3t_a = fse3t_a_crs(ji,jj,jk) 
     421                  !                                         ! tracer content at Before, now and after 
     422                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
     423                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n 
     424                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a 
     425                  ! 
     426                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b 
     427                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b 
     428                  ! 
     429                  ze3t_f = ze3t_n + atfp * ze3t_d 
     430                  ztc_f  = ztc_n  + atfp * ztc_d 
     431                  ! 
     432                  IF( jk == 1 ) THEN           ! first level  
     433                     ze3t_f = ze3t_f - zfact2 * ( emp_b_crs(ji,jj) - emp_crs(ji,jj) + rnf_crs(ji,jj) - rnf_b_crs(ji,jj) ) 
     434                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 
     435                  ENDIF 
     436!cbr as it is a subroutine dedicated to crs, TRA options are not necessary 
     437!cbr                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only) 
     438!cbr                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
     439!cbr 
     440!cbr                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )   &                  ! river runoffs 
     441!cbr                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
     442!cbr                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj) 
     443 
     444                  ze3t_f = 1.e0 / ze3t_f 
     445                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
     446                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
     447                  ! 
     448                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only) 
     449                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d ) 
     450                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
     451                  ENDIF 
     452               END DO 
     453            END DO 
     454         END DO 
     455         !  
     456      END DO 
     457      ! 
     458   END SUBROUTINE tra_nxt_vvl_crs 
     459 
     460 
    351461   !!====================================================================== 
    352462END MODULE tranxt 

  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp_crs.F90

    r5601 r6772  
    7878      !! ** Action  : - pta  becomes the after tracer 
    7979      !!--------------------------------------------------------------------- 
     80      USE ieee_arithmetic 
    8081      ! 
    8182      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
    … …  
    9091      REAL(wp) ::  zrhs, ze3tb, ze3tn, ze3ta   ! local scalars 
    9192      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt,zwd,zws 
     93      REAL(wp) ::  zmin,zmax 
    9294      !!--------------------------------------------------------------------- 
    9395      ! 
    … …  
    135137               END DO 
    136138            ELSE IF( l_traldf_rot ) THEN     ! standard isoneutral diff 
    137                DO jk = 2, jpkm1 
    138                   DO jj = 2, jpjm1 
    139                      DO ji = fs_2, fs_jpim1   ! vector opt. 
     139              DO jk = 2, jpkm1 
     140                  DO jj = 2, jpj_crs-1 
     141                     DO ji = 2, jpi_crs-1 
    140142                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk)                       & 
    141143                           &                          * (  wslpi_crs(ji,jj,jk) * wslpi_crs(ji,jj,jk)   & 
    … …  
    148150#endif 
    149151            DO jk = 1, jpkm1 
    150                DO jj = 2, jpjm1 
    151                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    152                      ze3ta =  ( 1. - r_vvl ) +        r_vvl   * e3t_crs(ji,jj,jk)   ! after scale factor at T-point 
    153                      ze3tn =         r_vvl   + ( 1. - r_vvl ) * e3t_crs(ji,jj,jk)   ! now   scale factor at T-point 
     152               DO jj = 2, jpj_crs-1 
     153                  DO ji = 2, jpi_crs-1 
     154 
     155#if defined key_vvl 
     156                     ze3ta =  ( 1. - r_vvl ) +        r_vvl   * fse3t_a_crs(ji,jj,jk)   ! after scale factor at T-point 
     157                     ze3tn =         r_vvl   + ( 1. - r_vvl ) * fse3t_n_crs(ji,jj,jk)   ! now   scale factor at T-point 
     158#else 
     159                     ze3ta =  ( 1. - r_vvl ) +        r_vvl   * e3t_0_crs(ji,jj,jk)   ! after scale factor at T-point 
     160                     ze3tn =         r_vvl   + ( 1. - r_vvl ) * e3t_0_crs(ji,jj,jk)   ! now   scale factor at T-point 
     161#endif 
    154162                     !cbr zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * e3w_1d(jk  ) )  !cc 
    155163                     !cbr zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * e3w_1d(jk+1) )  !cc 
    156                      zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * e3w_max_crs(ji,jj,jk) ) 
    157                      zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * e3w_max_crs(ji,jj,jk+1) ) 
     164                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w_max_crs(ji,jj,jk) ) 
     165                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w_max_crs(ji,jj,jk+1) ) 
    158166                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    159167                 END DO 
    … …  
    182190            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
    183191            ! done once for all passive tracers (so included in the IF instruction) 
    184             DO jj = 2, jpjm1 
    185                DO ji = fs_2, fs_jpim1 
     192            DO jj = 2, jpj_crs-1 
     193               DO ji = 2, jpi_crs-1 
    186194                  zwt(ji,jj,1) = zwd(ji,jj,1) 
    187195               END DO 
    188196            END DO 
    189197            DO jk = 2, jpkm1 
    190                DO jj = 2, jpjm1 
    191                   DO ji = fs_2, fs_jpim1 
     198               DO jj = 2, jpj_crs-1 
     199                  DO ji = 2, jpi_crs-1 
    192200                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
    193201                  END DO 
    … …  
    198206         !     
    199207         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    200          DO jj = 2, jpjm1 
    201             DO ji = fs_2, fs_jpim1 
    202                ze3tb = ( 1. - r_vvl ) + r_vvl * e3t_crs(ji,jj,1) 
    203                ze3tn = ( 1. - r_vvl ) + r_vvl * e3t_crs(ji,jj,1) 
     208         DO jj = 2, jpj_crs-1 
     209            DO ji = 2, jpi_crs-1 
     210#if defined key_vvl 
     211               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b_crs(ji,jj,1) 
     212               ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t_n_crs(ji,jj,1) 
     213#else 
     214               ze3tb = ( 1. - r_vvl ) + r_vvl * e3t_0_crs(ji,jj,1) 
     215               ze3tn = ( 1. - r_vvl ) + r_vvl * e3t_0_crs(ji,jj,1) 
     216#endif 
    204217               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
    205218            END DO 
    … …  
    207220 
    208221         DO jk = 2, jpkm1 
    209             DO jj = 2, jpjm1 
    210                DO ji = fs_2, fs_jpim1 
    211                   ze3tb = ( 1. - r_vvl ) + r_vvl * e3t_crs(ji,jj,jk) 
    212                   ze3tn = ( 1. - r_vvl ) + r_vvl * e3t_crs(ji,jj,jk) 
     222            DO jj = 2, jpj_crs-1 
     223               DO ji = 2, jpi_crs-1 
     224#if defined key_vvl 
     225                  ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b_crs(ji,jj,jk) 
     226                  ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t_n_crs(ji,jj,jk) 
     227#else 
     228                  ze3tb = ( 1. - r_vvl ) + r_vvl * e3t_0_crs(ji,jj,jk) 
     229                  ze3tn = ( 1. - r_vvl ) + r_vvl * e3t_0_crs(ji,jj,jk) 
     230#endif 
    213231                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side  
    214                   pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 
    215                  
     232                  pta(ji,jj,jk,jn) = zrhs  - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 
    216233               END DO 
    217234            END DO 
    … …  
    219236 
    220237         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
    221          DO jj = 2, jpjm1 
    222             DO ji = fs_2, fs_jpim1 
     238         DO jj = 2, jpj_crs-1 
     239            DO ji = 2, jpi_crs-1 
    223240               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask_crs(ji,jj,jpkm1) 
    224241            END DO 
    225242         END DO 
    226243         DO jk = jpk-2, 1, -1 
    227             DO jj = 2, jpjm1 
    228                DO ji = fs_2, fs_jpim1 
     244            DO jj = 2, jpj_crs-1 
     245               DO ji = 2, jpi_crs-1 
    229246                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
    230                      &             / zwt(ji,jj,jk) * tmask_crs(ji,jj,jk) 
    231                   
    232                END DO 
    233             END DO 
    234          END DO 
    235  
     247                    &             / zwt(ji,jj,jk) * tmask_crs(ji,jj,jk) 
     248              END DO 
     249            END DO 
     250         END DO 
    236251         !                                            ! ================= ! 
    237252      END DO                                          !  end tracer loop  ! 

  • branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde_crs.F90

    r5601 r6772  
    9696      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
    9797      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
    98   !cc    REAL(wp), POINTER, DIMENSION(:,:  ) ::  zri, zrj, zhi, zhj 
    99   !cc    REAL(wp), POINTER, DIMENSION(:,:,:) ::  zti, zte    ! interpolated value of tracer 
    10098      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::  zri, zrj, zhi, zhj 
    10199      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zti, zte    ! interpolated value of tracer 
    … …  
    105103      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_crs') 
    106104      ! 
    107 !!      CALL wrk_alloc( jpi, jpj,       zri, zrj, zhi, zhj )  
    108 !!      CALL wrk_alloc( jpi, jpj, kjpt, zti, zte           )  
    109105      ALLOCATE( zri(jpi_crs,jpj_crs) , zrj(jpi_crs,jpj_crs), zte(jpi_crs ,jpj_crs ,kjpt), & 
    110106         &      zhi(jpi_crs,jpj_crs) , zhj(jpi_crs,jpj_crs), zti(jpi_crs ,jpj_crs ,kjpt)) 
    … …  
    112108      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    113109         ! 
    114 # if defined key_vectopt_loop 
    115          jj = 1 
    116          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    117 # else 
    118110         DO jj = 1, jpjm1 
    119111            DO ji = 1, jpim1 
    120 # endif 
     112 
    121113               iku = mbku_crs(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    122114               ikv = mbkv_crs(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
    123           !     ze3wu = e3w_crs(ji+1,jj  ,iku) - e3w_crs(ji,jj,iku) 
    124           !     ze3wv = e3w_crs(ji  ,jj+1,ikv) - e3w_crs(ji,jj,ikv) 
    125                ze3wu = e3w_max_crs(ji+1,jj  ,iku) - e3w_max_crs(ji,jj,iku) 
    126                ze3wv = e3w_max_crs(ji  ,jj+1,ikv) - e3w_max_crs(ji,jj,ikv) 
     115               ze3wu = e3w_max_0_crs(ji+1,jj  ,iku) - e3w_max_0_crs(ji,jj,iku) 
     116               ze3wv = e3w_max_0_crs(ji  ,jj+1,ikv) - e3w_max_0_crs(ji,jj,ikv) 
    127117               ! 
    128118               ! i- direction 
    129119               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    130                   zmaxu =  ze3wu / e3w_max_crs(ji+1,jj,iku)   
    131                  !    zmaxu =  ze3wu / e3w_crs(ji+1,jj,iku) 
     120#if defined key_vvl 
     121                  zmaxu =  ze3wu / e3w_max_n_crs(ji+1,jj,iku)   
     122#else 
     123                  zmaxu =  ze3wu / e3w_max_0_crs(ji+1,jj,iku)   
     124#endif 
    132125                  ! interpolated values of tracers 
    133126                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
    … …  
    135128                  pgtu(ji,jj,jn) = umask_crs(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    136129               ELSE                           ! case 2 
    137                   zmaxu = -ze3wu / e3w_max_crs(ji,jj,iku) 
    138                  !    zmaxu = -ze3wu / e3w_crs(ji,jj,iku) 
     130#if defined key_vvl 
     131                  zmaxu = -ze3wu / e3w_max_n_crs(ji,jj,iku) 
     132#else 
     133                  zmaxu = -ze3wu / e3w_max_0_crs(ji,jj,iku) 
     134#endif 
    139135                  ! interpolated values of tracers 
    140136                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    … …  
    145141               ! j- direction 
    146142               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    147                   zmaxv =  ze3wv / e3w_max_crs(ji,jj+1,ikv) 
    148                !      zmaxv =  ze3wv / e3w_crs(ji,jj+1,ikv) 
     143#if defined key_vvl 
     144                  zmaxv =  ze3wv / e3w_max_n_crs(ji,jj+1,ikv) 
     145#else 
     146                  zmaxv =  ze3wv / e3w_max_0_crs(ji,jj+1,ikv) 
     147#endif 
    149148                  ! interpolated values of tracers 
    150149                  zte(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    … …  
    152151                  pgtv(ji,jj,jn) = vmask_crs(ji,jj,1) * ( zte(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    153152               ELSE                           ! case 2 
    154                   zmaxv =  -ze3wv / e3w_max_crs(ji,jj,ikv) 
    155                 !     zmaxv = -ze3wv / e3w_crs(ji,jj,ikv) 
     153#if defined key_vvl 
     154                  zmaxv =  -ze3wv / e3w_max_n_crs(ji,jj,ikv) 
     155#else 
     156                  zmaxv =  -ze3wv / e3w_max_0_crs(ji,jj,ikv) 
     157#endif 
    156158                  ! interpolated values of tracers 
    157159                  zte(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    … …  
    160162               ENDIF 
    161163 
    162 # if ! defined key_vectopt_loop 
    163164            END DO 
    164 # endif 
    165165         END DO 
    166166         CALL crs_lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL crs_lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    167167         ! 
    168168      END DO 
    169 !WRITE(numout,*) ' test24 ', e3w_max_crs 
     169 
    170170      ! horizontal derivative of density anomalies (rd) 
    171171      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    172 # if defined key_vectopt_loop 
    173          jj = 1 
    174          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    175 # else 
    176172         DO jj = 1, jpjm1 
    177173            DO ji = 1, jpim1 
    178 # endif 
     174 
    179175               iku = mbku_crs(ji,jj) 
    180176               ikv = mbkv_crs(ji,jj) 
    181    !cc             ze3wu  = e3w_max_crs(ji+1,jj  ,iku) - e3w_max_crs(ji,jj,iku)   !gradiant horizontal pas de max 
    182                ze3wu  = e3w_crs(ji+1,jj  ,iku) - e3w_crs(ji,jj,iku) 
    183        !cc        ze3wv  = e3w_max_crs(ji  ,jj+1,ikv) - e3w_max_crs(ji,jj,ikv) 
    184                ze3wv  = e3w_crs(ji  ,jj+1,ikv) - e3w_crs(ji,jj,ikv) 
    185                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_crs(ji  ,jj,iku)     ! i-direction: case 1 
    186                ELSE                        ;   zhi(ji,jj) = gdept_crs(ji+1,jj,iku)     ! -     -      case 2 
    187                ENDIF 
    188                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_crs(ji,jj  ,ikv)     ! j-direction: case 1 
    189                ELSE                        ;   zhj(ji,jj) = gdept_crs(ji,jj+1,ikv)     ! -     -      case 2 
    190                ENDIF 
    191 # if ! defined key_vectopt_loop 
     177               ze3wu  = e3w_0_crs(ji+1,jj  ,iku) - e3w_0_crs(ji,jj,iku) 
     178               ze3wv  = e3w_0_crs(ji  ,jj+1,ikv) - e3w_0_crs(ji,jj,ikv) 
     179               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept_crs(ji  ,jj,iku)     ! i-direction: case 1 
     180               ELSE                        ;   zhi(ji,jj) = fsdept_crs(ji+1,jj,iku)     ! -     -      case 2 
     181               ENDIF 
     182               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept_crs(ji,jj  ,ikv)     ! j-direction: case 1 
     183               ELSE                        ;   zhj(ji,jj) = fsdept_crs(ji,jj+1,ikv)     ! -     -      case 2 
     184               ENDIF 
     185 
    192186            END DO 
    193 # endif 
    194187         END DO 
    195188         CALL eos_crs( zti, zhi, zri )   
    196189         CALL eos_crs( zte, zhj, zrj ) 
     190 
    197191         ! Gradient of density at the last level  
    198 # if defined key_vectopt_loop 
    199          jj = 1 
    200          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    201 # else 
    202192         DO jj = 1, jpjm1 
    203193            DO ji = 1, jpim1 
    204 # endif 
    205194               iku = mbku_crs(ji,jj) 
    206195               ikv = mbkv_crs(ji,jj) 
    207       !         ze3wu  = e3w_max_crs(ji+1,jj  ,iku) - e3w_max_crs(ji,jj,iku)         gradient horizontal 
    208                 ze3wu  = e3w_crs(ji+1,jj  ,iku) - e3w_crs(ji,jj,iku) 
    209       !         ze3wv  = e3w_max_crs(ji  ,jj+1,ikv) - e3w_max_crs(ji,jj,ikv)         gradient horizontal 
    210                 ze3wv  = e3w_crs(ji  ,jj+1,ikv) - e3w_crs(ji,jj,ikv) 
     196                ze3wu  = e3w_0_crs(ji+1,jj  ,iku) - e3w_0_crs(ji,jj,iku) 
     197                ze3wv  = e3w_0_crs(ji  ,jj+1,ikv) - e3w_0_crs(ji,jj,ikv) 
    211198               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask_crs(ji,jj,1) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    212199               ELSE                        ;   pgru(ji,jj) = umask_crs(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
    … …  
    215202               ELSE                        ;   pgrv(ji,jj) = vmask_crs(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
    216203               ENDIF 
    217 # if ! defined key_vectopt_loop 
     204 
    218205            END DO 
    219 # endif 
    220206         END DO 
    221  
    222207 
    223208         CALL crs_lbc_lnk( pgru , 'U', -1. )   ;   CALL crs_lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
    … …  
    225210      END IF 
    226211      ! 
    227       !!ccCALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj )  
    228       !!ccCALL wrk_dealloc( jpi, jpj, kjpt, zti, zte           )  
    229212      DEALLOCATE( zri , zrj, zte, zhi, zhj, zti) 
    230213      ! 
Note: See TracChangeset for help on using the changeset viewer.