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 9528 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 – NEMO

Ignore:
Timestamp:
2018-04-30T15:39:18+02:00 (6 years ago)
Author:
gm
Message:

dev_merge_2017: add to new vorticity scheme (EET= EEN using e3t instead of e3f, and ENT= energy conserving at T-point). NB: not used in ref config.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r9506 r9528  
    3636   USE sbcapr          ! surface boundary condition: atmospheric pressure 
    3737   USE dynadv    , ONLY: ln_dynadv_vec 
     38   USE dynvor          ! vortivity scheme indicators 
    3839   USE phycst          ! physical constants 
    3940   USE dynvor          ! vorticity term 
     
    103104      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 
    104105      ! 
    105       IF( ln_dynvor_een )   ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
    106          &                            ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 
     106      IF( ln_dynvor_een .OR. ln_dynvor_eeT )   & 
     107         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
     108         &               ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 
    107109         ! 
    108110      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) ) 
     
    149151      INTEGER  ::   ikbu, iktu, noffset   ! local integers 
    150152      INTEGER  ::   ikbv, iktv            !   -      - 
    151       REAL(wp) ::   r1_2dt_b, z2dt_bf        ! local scalars 
    152       REAL(wp) ::   zx1, zx2, zu_spg, zhura  !   -      - 
    153       REAL(wp) ::   zy1, zy2, zv_spg, zhvra  !   -      - 
    154       REAL(wp) ::   za0, za1, za2, za3       !   -      - 
    155       REAL(wp) ::   zmdi, zztmp              !   -      - 
     153      REAL(wp) ::   r1_2dt_b, z2dt_bf               ! local scalars 
     154      REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      - 
     155      REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      - 
     156      REAL(wp) ::   za0, za1, za2, za3              !   -      - 
     157      REAL(wp) ::   zmdi, zztmp            , z1_ht  !   -      - 
    156158      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 
    157159      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 
    158160      REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 
    159       REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e 
     161      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 
    160162      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 
    161163      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points 
     
    237239      ! 
    238240      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 
    239          IF( ln_dynvor_een ) THEN               !==  EEN scheme  ==! 
     241         ! 
     242         SELECT CASE( nvor_scheme ) 
     243         CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme) 
    240244            SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
    241245            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     
    250254               DO jj = 1, jpjm1 
    251255                  DO ji = 1, jpim1 
    252                      zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     & 
    253                         &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   & 
    254                         &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
    255                         &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
     256!!gm bug for ocean cavities 
     257!                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     & 
     258!                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   & 
     259!                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
     260!                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
     261!! replace by : 
     262                     zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      & 
     263                        &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   & 
     264                        &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      & 
     265                        &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   ) 
     266!!gm end 
    256267                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    257268                  END DO 
     
    270281            END DO 
    271282            ! 
    272          ELSE                                !== all other schemes (ENE, ENS, MIX) 
     283         CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme) 
     284            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     285            DO jj = 2, jpj 
     286               DO ji = 2, jpi 
     287                  z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 
     288                  ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht 
     289                  ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht 
     290                  ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 
     291                  ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht 
     292               END DO 
     293            END DO 
     294            ! 
     295         CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT ! 
     296            ! 
    273297            zwz(:,:) = 0._wp 
    274298            zhf(:,:) = 0._wp 
     
    280304!!gm  
    281305!!             
    282               IF( .NOT.ln_sco ) THEN 
     306            IF( .NOT.ln_sco ) THEN 
    283307   
    284308   !!gm  agree the JC comment  : this should be done in a much clear way 
     
    290314   !              ENDIF 
    291315   !              zhf(:,:) = gdepw_0(:,:,jk+1) 
    292                ELSE 
    293                  !zhf(:,:) = hbatf(:,:) 
    294                  DO jj = 1, jpjm1 
    295                    DO ji = 1, jpim1 
    296                      zhf(ji,jj) = MAX( 0._wp,                                & 
    297                                 & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     & 
    298                                 &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     & 
    299                                 &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     & 
    300                                 &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   & 
    301                                 & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +& 
    302                                 &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +& 
    303                                 &   rsmall  ) ) 
    304                    END DO 
    305                  END DO 
    306               END IF 
    307    
    308               DO jj = 1, jpjm1 
    309                  zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
    310               END DO 
     316               ! 
     317            ELSE 
     318               ! 
     319               !zhf(:,:) = hbatf(:,:) 
     320               DO jj = 1, jpjm1 
     321                  DO ji = 1, jpim1 
     322!!gm Bug here in case of ocean cavities, further more ht_0 is masked by definition ==>> remove mask multiplication  
     323!                     zhf(ji,jj) = MAX( 0._wp,                                & 
     324!                                & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     & 
     325!                                &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     & 
     326!                                &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     & 
     327!                                &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   & 
     328!                                & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +& 
     329!                                &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +& 
     330!                                &   rsmall  ) ) 
     331!!gm replaced by : 
     332                     zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          & 
     333                        &              + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      & 
     334                        &       / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          & 
     335                        &              + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  ) 
    311336!!gm end 
    312  
     337                  END DO 
     338               END DO 
     339            ENDIF 
     340            ! 
     341            DO jj = 1, jpjm1 
     342               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     343            END DO 
     344            ! 
    313345            DO jk = 1, jpkm1 
    314346               DO jj = 1, jpjm1 
     
    324356            END DO 
    325357            zwz(:,:) = ff_f(:,:) * zwz(:,:) 
    326          ENDIF 
     358         END SELECT 
    327359      ENDIF 
    328360      ! 
     
    369401      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    370402      !                                   ! -------------------------------------------------------- 
     403      ! 
    371404      zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes  
    372405      zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 
    373406      ! 
    374       IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
     407      SELECT CASE( nvor_scheme ) 
     408      CASE( np_ENT )                ! enstrophy conserving scheme (f-point) 
     409         DO jj = 2, jpjm1 
     410            DO ji = 2, jpim1   ! vector opt. 
     411               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    & 
     412                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   & 
     413                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   ) 
     414                  ! 
     415               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    & 
     416                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   &  
     417                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   )  
     418            END DO   
     419         END DO   
     420         !          
     421      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX 
    375422         DO jj = 2, jpjm1 
    376423            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    385432         END DO 
    386433         ! 
    387       ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme 
     434      CASE( np_ENS )                ! enstrophy conserving scheme (f-point) 
    388435         DO jj = 2, jpjm1 
    389436            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    397444         END DO 
    398445         ! 
    399       ELSEIF ( ln_dynvor_een ) THEN  ! enstrophy and energy conserving scheme 
     446      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)          
    400447         DO jj = 2, jpjm1 
    401448            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    411458         END DO 
    412459         ! 
    413       ENDIF  
     460      END SELECT 
    414461      ! 
    415462      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    416463      !                                   ! ---------------------------------------------------- 
    417464      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    418         IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters 
    419            DO jj = 2, jpjm1 
    420               DO ji = 2, jpim1  
    421                 ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
    422                      &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
    423                      &    MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     465         IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters 
     466            DO jj = 2, jpjm1 
     467               DO ji = 2, jpim1  
     468                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
     469                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     470                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     471                     &                                                         > rn_wdmin1 + rn_wdmin2 
     472                  ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
     473                     &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     474                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     475                  IF(ll_tmp1) THEN 
     476                     zcpx(ji,jj) = 1.0_wp 
     477                  ELSEIF(ll_tmp2) THEN 
     478                     ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     479                     zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
     480                                 &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     481                     zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
     482                  ELSE 
     483                     zcpx(ji,jj) = 0._wp 
     484                  ENDIF 
     485                  ! 
     486                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                & 
     487                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
     488                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
    424489                     &                                                       > rn_wdmin1 + rn_wdmin2 
    425                 ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
    426                      &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
    427                      &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    428     
    429                 IF(ll_tmp1) THEN 
    430                   zcpx(ji,jj) = 1.0_wp 
    431                 ELSE IF(ll_tmp2) THEN 
    432                   ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    433                   zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    434                               &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
    435                   zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
    436  
    437                 ELSE 
    438                   zcpx(ji,jj) = 0._wp 
    439                 END IF 
    440           
    441                 ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                & 
    442                      &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
    443                      &    MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
    444                      &                                                       > rn_wdmin1 + rn_wdmin2 
    445                 ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
    446                      &    MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
    447                      &    MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     490                  ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
     491                     &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     492                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    448493   
    449                 IF(ll_tmp1) THEN 
    450                   zcpy(ji,jj) = 1.0_wp 
    451                 ELSE IF(ll_tmp2) THEN 
    452                   ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    453                   zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    454                               &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
    455                   zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 
    456  
    457                 ELSE 
    458                   zcpy(ji,jj) = 0._wp 
    459                 END IF 
    460               END DO 
    461            END DO 
    462   
    463            DO jj = 2, jpjm1 
    464               DO ji = 2, jpim1 
    465                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    466                         &                        * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
    467                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    468                         &                        * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth 
    469  
    470               END DO 
    471            END DO 
    472  
     494                  IF(ll_tmp1) THEN 
     495                     zcpy(ji,jj) = 1.0_wp 
     496                  ELSE IF(ll_tmp2) THEN 
     497                     ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     498                     zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
     499                        &             / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     500                     zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  ) 
     501                  ELSE 
     502                     zcpy(ji,jj) = 0._wp 
     503                  ENDIF 
     504               END DO 
     505            END DO 
     506            ! 
     507            DO jj = 2, jpjm1 
     508               DO ji = 2, jpim1 
     509                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     510                     &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
     511                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     512                     &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth 
     513               END DO 
     514            END DO 
     515            ! 
    473516         ELSE 
    474517            ! 
     
    675718      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
    676719      vn_adv(:,:) = 0._wp 
    677       !                                             ! ==================== ! 
    678  
    679       IF (ln_wd_dl) THEN 
     720      ! 
     721      IF( ln_wd_dl ) THEN 
    680722         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)  
    681723         zvwdmask(:,:) = 0._wp  !  
    682          zuwdav2(:,:) = 0._wp  
    683          zvwdav2(:,:) = 0._wp    
     724         zuwdav2 (:,:) = 0._wp  
     725         zvwdav2 (:,:) = 0._wp    
    684726      END IF  
    685727 
    686  
     728      !                                             ! ==================== ! 
    687729      DO jn = 1, icycle                             !  sub-time-step loop  ! 
    688730         !                                          ! ==================== ! 
     
    715757            ! set wetting & drying mask at tracer points for this barotropic sub-step  
    716758            IF ( ln_wd_dl ) THEN  
    717  
     759               ! 
    718760               IF ( ln_wd_dl_rmp ) THEN  
    719761                  DO jj = 1, jpj                                  
     
    736778                        ELSE  
    737779                           ztwdmask(ji,jj) = 0._wp 
    738                         END IF 
     780                        ENDIF 
    739781                     END DO 
    740782                  END DO 
    741                END IF  
    742  
    743             END IF  
    744             
    745  
     783               ENDIF  
     784               ! 
     785            ENDIF  
     786            ! 
    746787            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
    747788               DO ji = 2, fs_jpim1   ! Vector opt. 
     
    756797            CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 
    757798            ! 
    758             zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    759             zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     799            zhup2_e(:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
     800            zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 
     801            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 
    760802         ELSE 
    761             zhup2_e (:,:) = hu_n(:,:) 
    762             zhvp2_e (:,:) = hv_n(:,:) 
     803            zhup2_e(:,:) = hu_n(:,:) 
     804            zhvp2_e(:,:) = hv_n(:,:) 
     805            zhtp2_e(:,:) = ht_n(:,:) 
    763806         ENDIF 
    764807         !                                                !* after ssh 
     
    795838         ENDIF 
    796839#endif 
    797          IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
     840         IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    798841 
    799842         IF ( ln_wd_dl ) THEN  
    800  
    801  
    802 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells  
    803  
     843            ! 
     844            ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells  
     845            ! 
    804846            DO jj = 1, jpjm1                                  
    805847               DO ji = 1, jpim1    
     
    821863               END DO 
    822864            END DO 
    823  
    824  
    825          END IF     
     865            ! 
     866         ENDIF     
    826867          
    827868         ! Sum over sub-time-steps to compute advective velocities 
     
    896937         ENDIF 
    897938         ! 
    898          zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
    899           &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
     939         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   & 
     940            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
     941          
    900942         IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters 
    901943           DO jj = 2, jpjm1 
     
    917959                ELSE 
    918960                  zcpx(ji,jj) = 0._wp 
    919                 END IF 
    920           
     961                ENDIF 
     962                ! 
    921963                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
    922964                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            & 
     
    929971                IF(ll_tmp1) THEN 
    930972                  zcpy(ji,jj) = 1.0_wp 
    931                 ELSE IF(ll_tmp2) THEN 
     973                ELSEIF(ll_tmp2) THEN 
    932974                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here 
    933975                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 
     
    935977                ELSE 
    936978                  zcpy(ji,jj) = 0._wp 
    937                 END IF 
     979                ENDIF 
    938980              END DO 
    939981           END DO 
    940          END IF 
     982         ENDIF 
    941983         ! 
    942984         ! Compute associated depths at U and V points: 
     
    955997               END DO 
    956998            END DO 
    957  
     999            ! 
    9581000         ENDIF 
    9591001         ! 
     
    9651007         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    9661008         ! 
    967          IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN     !==  energy conserving or mixed scheme  ==! 
     1009         SELECT CASE( nvor_scheme ) 
     1010         CASE( np_ENT )             ! energy conserving scheme (t-point) 
     1011         DO jj = 2, jpjm1 
     1012            DO ji = 2, jpim1   ! vector opt. 
     1013 
     1014               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 
     1015               z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
     1016             
     1017               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   & 
     1018                  &               * (  e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) )   & 
     1019                  &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   ) 
     1020                  ! 
     1021               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    & 
     1022                  &               * (  e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) )   &  
     1023                  &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   )  
     1024            END DO   
     1025         END DO   
     1026         !          
     1027         CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point) 
    9681028            DO jj = 2, jpjm1 
    9691029               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    9771037            END DO 
    9781038            ! 
    979          ELSEIF ( ln_dynvor_ens ) THEN                   !==  enstrophy conserving scheme  ==! 
     1039         CASE( np_ENS )             ! enstrophy conserving scheme (f-point) 
    9801040            DO jj = 2, jpjm1 
    9811041               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    9891049            END DO 
    9901050            ! 
    991          ELSEIF ( ln_dynvor_een ) THEN                   !==  energy and enstrophy conserving scheme  ==! 
     1051         CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f) 
    9921052            DO jj = 2, jpjm1 
    9931053               DO ji = fs_2, fs_jpim1   ! vector opt. 
    994                   zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
    995                      &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
    996                      &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1) &  
    997                      &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    998                   zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &  
    999                      &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
    1000                      &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &  
    1001                      &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     1054                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  & 
     1055                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  & 
     1056                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  &  
     1057                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  ) 
     1058                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  &  
     1059                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  & 
     1060                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  &  
     1061                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  ) 
    10021062               END DO 
    10031063            END DO 
    10041064            !  
    1005          ENDIF 
     1065         END SELECT 
    10061066         ! 
    10071067         ! Add tidal astronomical forcing if defined 
     
    13411401      !! ** Purpose : Read or write time-splitting arrays in restart file 
    13421402      !!---------------------------------------------------------------------- 
    1343       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    1344       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    1345       ! 
    1346       INTEGER ::   id1, id2   ! local integers 
     1403      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     1404      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    13471405      !!---------------------------------------------------------------------- 
    13481406      ! 
Note: See TracChangeset for help on using the changeset viewer.