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 6808 for branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 – NEMO

Ignore:
Timestamp:
2016-07-19T10:38:35+02:00 (8 years ago)
Author:
jamesharle
Message:

merge with trunk@6232 for consistency with SSB code

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r5467 r6808  
    2828   USE sbc_oce         ! surface boundary condition: ocean 
    2929   USE sbcrnf          ! river runoffs 
     30   USE sbcisf          ! ice shelf melting 
    3031   USE zdf_oce         ! ocean vertical mixing 
    3132   USE domvvl          ! variable volume 
    32    USE dynspg_oce      ! surface     pressure gradient variables 
    33    USE dynhpg          ! hydrostatic pressure gradient  
    3433   USE trd_oce         ! trends: ocean variables 
    3534   USE trdtra          ! trends manager: tracers  
    3635   USE traqsr          ! penetrative solar radiation (needed for nksr) 
    3736   USE phycst          ! physical constant 
    38    USE ldftra_oce      ! lateral physics on tracers 
     37   USE ldftra          ! lateral physics on tracers 
     38   USE ldfslp 
    3939   USE bdy_oce         ! BDY open boundary condition variables 
    4040   USE bdytra          ! open boundary condition (bdy_tra routine) 
     
    4646   USE timing          ! Timing 
    4747#if defined key_agrif 
    48    USE agrif_opa_update 
    4948   USE agrif_opa_interp 
    5049#endif 
     
    5756   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt 
    5857 
    59    REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg 
    60  
    6158   !! * Substitutions 
    62 #  include "domzgr_substitute.h90" 
     59#  include "vectopt_loop_substitute.h90" 
    6360   !!---------------------------------------------------------------------- 
    6461   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)  
     
    8885      !!             domains (lk_agrif=T) 
    8986      !! 
    90       !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
    91       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     87      !! ** Action  : - tsb & tsn ready for the next time step 
    9288      !!---------------------------------------------------------------------- 
    9389      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    9490      !! 
    95       INTEGER  ::   jk, jn    ! dummy loop indices 
    96       REAL(wp) ::   zfact     ! local scalars 
     91      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     92      REAL(wp) ::   zfact            ! local scalars 
    9793      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    9894      !!---------------------------------------------------------------------- 
     
    104100         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 
    105101         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    106          ! 
    107          rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp)      ! Brown & Campana parameter for semi-implicit hpg 
    108102      ENDIF 
    109103 
    110104      ! Update after tracer on domain lateral boundaries 
    111105      !  
     106#if defined key_agrif 
     107      CALL Agrif_tra                     ! AGRIF zoom boundaries 
     108#endif 
     109      ! 
    112110      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign) 
    113111      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp ) 
     
    116114      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries 
    117115#endif 
    118 #if defined key_agrif 
    119       CALL Agrif_tra                     ! AGRIF zoom boundaries 
    120 #endif 
    121116  
    122117      ! set time step size (Euler/Leapfrog) 
    123       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler) 
    124       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog) 
     118      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =     rdt      ! at nit000             (Euler) 
     119      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog) 
    125120      ENDIF 
    126121 
     
    142137            END DO 
    143138         END DO 
     139         ! 
    144140      ELSE                                            ! Leap-Frog + Asselin filter time stepping 
    145141         ! 
    146          IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   & 
    147            &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)  
    148          ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level  
     142         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface  
     143         ELSE                   ;   CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa,   & 
     144           &                                                                sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface 
    149145         ENDIF 
    150       ENDIF  
    151       ! 
    152 #if defined key_agrif 
    153       ! Update tracer at AGRIF zoom boundaries 
    154       IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only 
    155 #endif       
    156       ! 
    157       ! trends computation 
     146         ! 
     147         DO jn = 1, jpts 
     148            CALL lbc_lnk( tsb(:,:,:,jn), 'T', 1._wp )  
     149            CALL lbc_lnk( tsn(:,:,:,jn), 'T', 1._wp ) 
     150            CALL lbc_lnk( tsa(:,:,:,jn), 'T', 1._wp ) 
     151         END DO 
     152      ENDIF      
     153      ! 
    158154      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt      
    159155         DO jk = 1, jpkm1 
    160             zfact = 1._wp / r2dtra(jk)              
     156            zfact = 1._wp / r2dt              
    161157            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact 
    162158            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact 
     
    184180      !!  
    185181      !! ** Method  : - Apply a Asselin time filter on now fields. 
    186       !!              - save in (ta,sa) an average over the three time levels  
    187       !!             which will be used to compute rdn and thus the semi-implicit 
    188       !!             hydrostatic pressure gradient (ln_dynhpg_imp = T) 
    189182      !!              - swap tracer fields to prepare the next time_step. 
    190       !!                This can be summurized for tempearture as: 
    191       !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T 
    192       !!             ztm = 0                                   otherwise 
    193       !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp) 
    194       !!             tb  = tn + atfp*[ tb - 2 tn + ta ] 
    195       !!             tn  = ta   
    196       !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call) 
    197       !! 
    198       !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
    199       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    200       !!---------------------------------------------------------------------- 
    201       INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index 
    202       INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index 
    203       CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator) 
    204       INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers 
    205       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields 
    206       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields 
    207       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend 
     183      !! 
     184      !! ** Action  : - tsb & tsn ready for the next time step 
     185      !!---------------------------------------------------------------------- 
     186      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index 
     187      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index 
     188      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator) 
     189      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers 
     190      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields 
     191      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields 
     192      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend 
    208193      ! 
    209194      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    210       LOGICAL  ::   ll_tra_hpg       ! local logical 
    211195      REAL(wp) ::   ztn, ztd         ! local scalars 
    212196      !!---------------------------------------------------------------------- 
    213  
     197      ! 
    214198      IF( kt == kit000 )  THEN 
    215199         IF(lwp) WRITE(numout,*) 
     
    218202      ENDIF 
    219203      ! 
    220       IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg     
    221       ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg 
    222       ENDIF 
    223       ! 
    224204      DO jn = 1, kjpt 
    225205         ! 
    226206         DO jk = 1, jpkm1 
    227             DO jj = 1, jpj 
    228                DO ji = 1, jpi 
     207            DO jj = 2, jpjm1 
     208               DO ji = fs_2, fs_jpim1 
    229209                  ztn = ptn(ji,jj,jk,jn)                                     
    230                   ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers 
    231                   ! 
    232                   ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn  
    233                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta 
    234                   ! 
    235                   IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average 
     210                  ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn)  ! time laplacian on tracers 
     211                  ! 
     212                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                      ! ptb <-- filtered ptn  
     213                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                      ! ptn <-- pta 
    236214               END DO 
    237215           END DO 
     
    251229      !!  
    252230      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
    253       !!              - save in (ta,sa) a thickness weighted average over the three  
    254       !!             time levels which will be used to compute rdn and thus the semi- 
    255       !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 
    256231      !!              - swap tracer fields to prepare the next time_step. 
    257       !!                This can be summurized for tempearture as: 
    258       !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T 
    259       !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )    
    260       !!             ztm = 0                                                       otherwise 
    261232      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
    262233      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] ) 
    263234      !!             tn  = ta  
    264       !!             ta  = zt        (NB: reset to 0 after eos_bn2 call) 
    265       !! 
    266       !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
    267       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    268       !!---------------------------------------------------------------------- 
    269       INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index 
    270       INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index 
    271       REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step 
    272       CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator) 
    273       INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers 
    274       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields 
    275       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields 
    276       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend 
    277       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content 
    278       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content 
    279  
    280       !!      
    281       LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical 
     235      !! 
     236      !! ** Action  : - tsb & tsn ready for the next time step 
     237      !!---------------------------------------------------------------------- 
     238      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index 
     239      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index 
     240      REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step 
     241      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator) 
     242      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers 
     243      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields 
     244      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields 
     245      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend 
     246      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc   ! surface tracer content 
     247      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc_b ! before surface tracer content 
     248      ! 
     249      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical 
    282250      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
    283251      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar 
     
    292260      ! 
    293261      IF( cdtype == 'TRA' )  THEN    
    294          ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg 
    295262         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration 
    296263         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs 
    297       ELSE                           
    298          ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg 
    299          ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration 
    300          ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs 
     264         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting 
     265      ELSE                          ! passive tracers case 
     266         ll_traqsr  = .FALSE.          ! NO solar penetration 
     267         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?   
     268         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??  
    301269      ENDIF 
    302270      ! 
    303271      DO jn = 1, kjpt       
    304272         DO jk = 1, jpkm1 
    305             zfact1 = atfp * p2dt(jk) 
    306             zfact2 = zfact1 / rau0 
    307             DO jj = 1, jpj 
    308                DO ji = 1, jpi 
    309                   ze3t_b = fse3t_b(ji,jj,jk) 
    310                   ze3t_n = fse3t_n(ji,jj,jk) 
    311                   ze3t_a = fse3t_a(ji,jj,jk) 
     273            zfact1 = atfp * p2dt 
     274            zfact2 = zfact1 * r1_rau0 
     275            DO jj = 2, jpjm1 
     276               DO ji = fs_2, fs_jpim1 
     277                  ze3t_b = e3t_b(ji,jj,jk) 
     278                  ze3t_n = e3t_n(ji,jj,jk) 
     279                  ze3t_a = e3t_a(ji,jj,jk) 
    312280                  !                                         ! tracer content at Before, now and after 
    313281                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
     
    321289                  ztc_f  = ztc_n  + atfp * ztc_d 
    322290                  ! 
    323                   IF( jk == 1 ) THEN           ! first level  
    324                      ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) ) 
     291                  IF( jk == mikt(ji,jj) ) THEN           ! first level  
     292                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  & 
     293                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  & 
     294                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  ) 
    325295                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 
    326296                  ENDIF 
    327  
    328                   IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only) 
     297                  ! 
     298                  ! solar penetration (temperature only) 
     299                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            &  
    329300                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )  
    330  
    331                   IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )   &            ! river runoffs 
     301                     ! 
     302                  ! river runoff 
     303                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          & 
    332304                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) &  
    333                      &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj) 
    334  
     305                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj) 
     306                     ! 
     307                  ! ice shelf 
     308                  IF( ll_isf ) THEN 
     309                     ! level fully include in the Losch_2008 ice shelf boundary layer 
     310                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          & 
     311                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  & 
     312                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) 
     313                     ! level partially include in Losch_2008 ice shelf boundary layer  
     314                     IF ( jk == misfkb(ji,jj) )                                                   & 
     315                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  & 
     316                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj) 
     317                  END IF 
     318                  ! 
    335319                  ze3t_f = 1.e0 / ze3t_f 
    336320                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
    337321                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
    338322                  ! 
    339                   IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only) 
    340                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d ) 
    341                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
    342                   ENDIF 
    343323               END DO 
    344324            END DO 
Note: See TracChangeset for help on using the changeset viewer.