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

Changeset 3008


Ignore:
Timestamp:
2011-10-27T13:10:32+02:00 (12 years ago)
Author:
acc
Message:

Branch dev_NOC_2011_MERGE. #874. Step 6: Merge in changes from 2011/dev_r2802_NOCS_vvlfix branch

Location:
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r2779 r3008  
    2424   PRIVATE 
    2525 
    26    PUBLIC   dom_vvl       ! called by domain.F90 
    27    PUBLIC   dom_vvl_alloc ! called by nemogcm.F90 
    28  
    29    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ee_t, ee_u, ee_v, ee_f   !: ??? 
    30    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: ???  
     26   PUBLIC   dom_vvl         ! called by domain.F90 
     27   PUBLIC   dom_vvl_2       ! called by domain.F90 
     28   PUBLIC   dom_vvl_alloc   ! called by nemogcm.F90 
     29 
     30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: 1/H_0 at t-,u-,v-,f-points  
    3131 
    3232   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra  
    … …  
    4949      ! 
    5050      ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) ,     & 
    51          &      ee_t(jpi,jpj)     , ee_u(jpi,jpj)     , ee_v(jpi,jpj)     , ee_f(jpi,jpj)     ,     & 
    5251         &      r2dt        (jpk)                                                             , STAT=dom_vvl_alloc ) 
    5352         ! 
    … …  
    6261      !!                ***  ROUTINE dom_vvl  *** 
    6362      !!                    
    64       !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread 
    65       !!               ssh over the whole water column (scale factors) 
     63      !! ** Purpose :   compute mu coefficients at t-, u-, v- and f-points to  
     64      !!              spread ssh over the whole water column (scale factors) 
     65      !!                set the before and now ssh at u- and v-points  
     66      !!              (also f-point in now case) 
    6667      !!---------------------------------------------------------------------- 
    6768      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    68       USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3     ! 2D workspace 
     69      USE wrk_nemo, ONLY:   zee_t => wrk_2d_1, zee_u => wrk_2d_2, zee_v => wrk_2d_3, zee_f => wrk_2d_4   ! 2D workspace 
    6970      ! 
    7071      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    71       REAL(wp) ::   zcoefu , zcoefv   , zcoeff                   ! local scalars 
    72       REAL(wp) ::   zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1   !   -      - 
    73       !!---------------------------------------------------------------------- 
    74  
    75       IF( wrk_in_use(2, 1,2,3) ) THEN 
     72      REAL(wp) ::   zcoefu, zcoefv , zcoeff                ! local scalars 
     73      REAL(wp) ::   zvt   , zvt_ip1, zvt_jp1, zvt_ip1jp1   !   -      - 
     74      !!---------------------------------------------------------------------- 
     75 
     76      IF( wrk_in_use(2, 1,2,3,4) ) THEN 
    7677         CALL ctl_stop('dom_vvl: requested workspace arrays unavailable')   ;   RETURN 
    7778      ENDIF 
    … …  
    9798 
    9899      !                                 !==  mu computation  ==! 
    99       ee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level 
    100       ee_u(:,:) = fse3u_0(:,:,1) 
    101       ee_v(:,:) = fse3v_0(:,:,1) 
    102       ee_f(:,:) = fse3f_0(:,:,1) 
     100      zee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level 
     101      zee_u(:,:) = fse3u_0(:,:,1) 
     102      zee_v(:,:) = fse3v_0(:,:,1) 
     103      zee_f(:,:) = fse3f_0(:,:,1) 
    103104      DO jk = 2, jpkm1                          ! Sum of the masked vertical scale factors 
    104          ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 
    105          ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
    106          ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
     105         zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 
     106         zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
     107         zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
    107108         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask 
    108             ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk) 
     109            zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk) 
    109110         END DO 
    110111      END DO   
    111112      !                                         ! Compute and mask the inverse of the local depth at T, U, V and F points 
    112       ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1) 
    113       ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1) 
    114       ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1) 
     113      zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1) 
     114      zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1) 
     115      zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1) 
    115116      DO jj = 1, jpjm1                               ! f-point case fmask cannot be used  
    116          ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 
    117       END DO 
    118       CALL lbc_lnk( ee_f, 'F', 1. )                  ! lateral boundary condition on ee_f 
     117         zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 
     118      END DO 
     119      CALL lbc_lnk( zee_f, 'F', 1. )                 ! lateral boundary condition on ee_f 
    119120      ! 
    120121      DO jk = 1, jpk                            ! mu coefficients 
    121          mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels 
    122          muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk)     ! U-point at T levels 
    123          muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels 
     122         mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels 
     123         muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk)     ! U-point at T levels 
     124         muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels 
    124125      END DO 
    125126      DO jk = 1, jpk                                 ! F-point : fmask=shlat at coasts, use the product of umask 
    126127         DO jj = 1, jpjm1 
    127                muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels 
    128          END DO 
    129          muf(:,jpj,jk) = 0.e0 
     128               muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels 
     129         END DO 
     130         muf(:,jpj,jk) = 0._wp 
    130131      END DO 
    131132      CALL lbc_lnk( muf, 'F', 1. )                   ! lateral boundary condition 
    … …  
    139140      END DO 
    140141       
    141       ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 
    142       ! for ssh and scale factors 
    143       zs_t  (:,:) =         e1t(:,:) * e2t(:,:) 
    144       zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) ) 
    145       zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) ) 
    146  
    147142      DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points 
    148143         DO ji = 1, jpim1   ! NO vector opt. 
    149             zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj) 
    150             zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj) 
    151             zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) ) 
    152             ! before fields 
    153             zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  ) 
    154             zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  ) 
    155             zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1) 
    156             sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 
    157             sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 
    158             ! now fields 
    159             zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  ) 
    160             zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  ) 
    161             zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1) 
    162             zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1) 
    163             sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 
    164             sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 
    165             sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 ) 
     144            zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1) 
     145            zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1) 
     146            zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1) 
     147            ! 
     148            zvt           = e1e2t(ji  ,jj  ) * sshb(ji  ,jj  )    ! before fields 
     149            zvt_ip1       = e1e2t(ji+1,jj  ) * sshb(ji+1,jj  ) 
     150            zvt_jp1       = e1e2t(ji  ,jj+1) * sshb(ji  ,jj+1) 
     151            sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 
     152            sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 
     153            ! 
     154            zvt           = e1e2t(ji  ,jj  ) * sshn(ji  ,jj  )    ! now fields 
     155            zvt_ip1       = e1e2t(ji+1,jj  ) * sshn(ji+1,jj  ) 
     156            zvt_jp1       = e1e2t(ji  ,jj+1) * sshn(ji  ,jj+1) 
     157            zvt_ip1jp1    = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) 
     158            sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 
     159            sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 
     160            sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 ) 
    166161         END DO 
    167162      END DO 
    … …  
    169164      CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. ) 
    170165      CALL lbc_lnk( sshf_n, 'F', 1. ) 
    171  
    172                                                 ! initialise before scale factors at (u/v)-points 
    173       ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    174       DO jk = 1, jpkm1 
    175          DO jj = 1, jpjm1 
    176             DO ji = 1, jpim1 
    177                zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    178                zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    179                zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    180                fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    181                fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
    182             END DO 
    183          END DO 
    184       END DO 
    185       CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions 
    186       CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
    187       ! Add initial scale factor to scale factor anomaly 
    188       fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
    189       fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 
    190       ! 
    191       IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays') 
     166      ! 
     167      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays') 
    192168      ! 
    193169   END SUBROUTINE dom_vvl 
    194170 
     171 
     172   SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b ) 
     173      !!---------------------------------------------------------------------- 
     174      !!                ***  ROUTINE dom_vvl_2  *** 
     175      !!                    
     176      !! ** Purpose :   compute the vertical scale factors at u- and v-points 
     177      !!              in variable volume case. 
     178      !! 
     179      !! ** Method  :   In variable volume case (non linear sea surface) the  
     180      !!              the vertical scale factor at velocity points is computed 
     181      !!              as the average of the cell surface weighted e3t. 
     182      !!                It uses the sea surface heigth so it have to be initialized 
     183      !!              after ssh is read/set 
     184      !!---------------------------------------------------------------------- 
     185      INTEGER                   , INTENT(in   ) ::   kt               ! ocean time-step index 
     186      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3u_b, pe3v_b   ! before vertical scale factor at u- & v-pts 
     187      ! 
     188      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     189      INTEGER  ::   iku, ikv     ! local integers     
     190      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers 
     191      REAL(wp) ::   zvt          ! local scalars 
     192      !!---------------------------------------------------------------------- 
     193 
     194      IF( lwp .AND. kt == nit000 ) THEN 
     195         WRITE(numout,*) 
     196         WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization' 
     197         WRITE(numout,*) '~~~~~~~~~ ' 
     198         pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk) 
     199         pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk) 
     200      ENDIF 
     201       
     202      DO jk = 1, jpkm1           ! set the before scale factors at u- & v-points 
     203         DO jj = 2, jpjm1 
     204            DO ji = fs_2, fs_jpim1 
     205               zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj) 
     206               pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     207               pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     208            END DO 
     209         END DO 
     210      END DO 
     211 
     212      ! Correct scale factors at locations that have been individually modified in domhgr 
     213      ! Such modifications break the relationship between e1e2t and e1u*e2u etc. Recompute 
     214      ! scale factors ignoring the modified metric. 
     215      !                                                ! ===================== 
     216      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
     217         !                                             ! ===================== 
     218         IF( nn_cla == 0 ) THEN 
     219            ! 
     220            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified) 
     221            ij0 = 102   ;   ij1 = 102    
     222            DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     223               DO jj = mj0(ij0), mj1(ij1) 
     224                  DO ji = mi0(ii0), mi1(ii1) 
     225                     zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     226                     pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     227                  END DO 
     228               END DO 
     229            END DO 
     230            ! 
     231            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified) 
     232            ij0 =  88   ;   ij1 =  88    
     233            DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     234               DO jj = mj0(ij0), mj1(ij1) 
     235                  DO ji = mi0(ii0), mi1(ii1) 
     236                     zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     237                     pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     238                  END DO 
     239               END DO 
     240            END DO 
     241            DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     242               DO jj = mj0(ij0), mj1(ij1) 
     243                  DO ji = mi0(ii0), mi1(ii1) 
     244                     zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     245                     pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     246                  END DO 
     247               END DO 
     248            END DO 
     249         ENDIF 
     250 
     251         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified) 
     252         ij0 = 116   ;   ij1 = 116    
     253         DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     254            DO jj = mj0(ij0), mj1(ij1) 
     255               DO ji = mi0(ii0), mi1(ii1) 
     256                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     257                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     258               END DO 
     259            END DO 
     260         END DO 
     261         ! 
     262      ENDIF 
     263         !                                             ! ===================== 
     264      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
     265         !                                             ! ===================== 
     266 
     267         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified) 
     268         ij0 = 200   ;   ij1 = 200    
     269         DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     270            DO jj = mj0(ij0), mj1(ij1) 
     271               DO ji = mi0(ii0), mi1(ii1) 
     272                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     273                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     274               END DO 
     275            END DO 
     276         END DO 
     277 
     278         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified) 
     279         ij0 = 208   ;   ij1 = 208    
     280         DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     281            DO jj = mj0(ij0), mj1(ij1) 
     282               DO ji = mi0(ii0), mi1(ii1) 
     283                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     284                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     285               END DO 
     286            END DO 
     287         END DO 
     288 
     289         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified) 
     290         ij0 = 124   ;   ij1 = 125    
     291         DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     292            DO jj = mj0(ij0), mj1(ij1) 
     293               DO ji = mi0(ii0), mi1(ii1) 
     294                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     295                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     296               END DO 
     297            END DO 
     298         END DO 
     299 
     300         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 
     301         ij0 = 124   ;   ij1 = 125    
     302         DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     303            DO jj = mj0(ij0), mj1(ij1) 
     304               DO ji = mi0(ii0), mi1(ii1) 
     305                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     306                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     307               END DO 
     308            END DO 
     309         END DO 
     310 
     311         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified) 
     312         ij0 = 124   ;   ij1 = 125    
     313         DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     314            DO jj = mj0(ij0), mj1(ij1) 
     315               DO ji = mi0(ii0), mi1(ii1) 
     316                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     317                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     318               END DO 
     319            END DO 
     320         END DO 
     321 
     322         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified) 
     323         ij0 = 124   ;   ij1 = 125    
     324         DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     325            DO jj = mj0(ij0), mj1(ij1) 
     326               DO ji = mi0(ii0), mi1(ii1) 
     327                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     328                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     329               END DO 
     330            END DO 
     331         END DO 
     332 
     333         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified) 
     334         ij0 = 141   ;   ij1 = 142    
     335         DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     336            DO jj = mj0(ij0), mj1(ij1) 
     337               DO ji = mi0(ii0), mi1(ii1) 
     338                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     339                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     340               END DO 
     341            END DO 
     342         END DO 
     343 
     344         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified) 
     345         ij0 = 141   ;   ij1 = 142    
     346         DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     347            DO jj = mj0(ij0), mj1(ij1) 
     348               DO ji = mi0(ii0), mi1(ii1) 
     349                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     350                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     351               END DO 
     352            END DO 
     353         END DO 
     354 
     355         ! 
     356      ENDIF 
     357      !                                                ! ====================== 
     358      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration 
     359         !                                             ! ====================== 
     360         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified) 
     361         ij0 = 327   ;   ij1 = 327    
     362         DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     363            DO jj = mj0(ij0), mj1(ij1) 
     364               DO ji = mi0(ii0), mi1(ii1) 
     365                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     366                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     367               END DO 
     368            END DO 
     369         END DO 
     370         ! 
     371         ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u was modified) 
     372         ij0 = 343   ;   ij1 = 343    
     373         DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     374            DO jj = mj0(ij0), mj1(ij1) 
     375               DO ji = mi0(ii0), mi1(ii1) 
     376                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     377                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     378               END DO 
     379            END DO 
     380         END DO 
     381         ! 
     382         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified) 
     383         ij0 = 232   ;   ij1 = 232    
     384         DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     385            DO jj = mj0(ij0), mj1(ij1) 
     386               DO ji = mi0(ii0), mi1(ii1) 
     387                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     388                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     389               END DO 
     390            END DO 
     391         END DO 
     392         ! 
     393         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified) 
     394         ij0 = 232   ;   ij1 = 232    
     395         DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     396            DO jj = mj0(ij0), mj1(ij1) 
     397               DO ji = mi0(ii0), mi1(ii1) 
     398                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     399                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     400               END DO 
     401            END DO 
     402         END DO 
     403         ! 
     404         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified) 
     405         ij0 = 270   ;   ij1 = 270    
     406         DO jk = 1, jpkm1                 ! set the before scale factors at u-points 
     407            DO jj = mj0(ij0), mj1(ij1) 
     408               DO ji = mi0(ii0), mi1(ii1) 
     409                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 
     410                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 
     411               END DO 
     412            END DO 
     413         END DO 
     414         ! 
     415         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified) 
     416         ij0 = 232   ;   ij1 = 233    
     417         DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     418            DO jj = mj0(ij0), mj1(ij1) 
     419               DO ji = mi0(ii0), mi1(ii1) 
     420                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     421                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     422               END DO 
     423            END DO 
     424         END DO 
     425         ! 
     426         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified) 
     427         ij0 = 276   ;   ij1 = 276    
     428         DO jk = 1, jpkm1                 ! set the before scale factors at v-points 
     429            DO jj = mj0(ij0), mj1(ij1) 
     430               DO ji = mi0(ii0), mi1(ii1) 
     431                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 
     432                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 
     433               END DO 
     434            END DO 
     435         END DO 
     436         ! 
     437      ENDIF 
     438      ! End of individual corrections to scale factors 
     439 
     440      IF( ln_zps ) THEN          ! minimum of the e3t at partial cell level 
     441         DO jj = 2, jpjm1 
     442            DO ji = fs_2, fs_jpim1 
     443               iku = mbku(ji,jj) 
     444               ikv = mbkv(ji,jj) 
     445               pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj  ,iku) )  
     446               pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji  ,jj+1,ikv) )  
     447            END DO 
     448         END DO 
     449      ENDIF 
     450 
     451      pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:)      ! anomaly to avoid zero along closed boundary/extra halos 
     452      pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:) 
     453      CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions 
     454      CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. ) 
     455      pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:)      ! recover the full scale factor 
     456      pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:) 
     457      ! 
     458   END SUBROUTINE dom_vvl_2 
     459    
    195460#else 
    196461   !!---------------------------------------------------------------------- 
    … …  
    200465   SUBROUTINE dom_vvl 
    201466   END SUBROUTINE dom_vvl 
     467   SUBROUTINE dom_vvl_2(kdum, pudum, pvdum ) 
     468      USE par_kind 
     469      INTEGER                   , INTENT(in   ) ::   kdum        
     470      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pudum, pvdum 
     471   END SUBROUTINE dom_vvl_2 
    202472#endif 
    203473 

  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r2777 r3008  
    8383         neuler = 1                              ! Set time-step indicator at nit000 (leap-frog) 
    8484         CALL rst_read                           ! Read the restart file 
     85         !                                       ! define e3u_b, e3v_b from e3t_b read in restart file 
     86         CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 
    8587         CALL tra_swap                           ! swap 3D arrays (t,s)  in a 4D array (ts) 
    8688         CALL day_init                           ! model calendar (using both namelist and restart infos) 
    … …  
    9294         CALL day_init                           ! model calendar (using both namelist and restart infos) 
    9395         !                                       ! Initialization of ocean to zero 
    94          !   before fields     !       now fields      
    95          sshb (:,:)   = 0.e0   ;   sshn (:,:)   = 0.e0 
    96          ub   (:,:,:) = 0.e0   ;   un   (:,:,:) = 0.e0 
    97          vb   (:,:,:) = 0.e0   ;   vn   (:,:,:) = 0.e0   
    98          rotb (:,:,:) = 0.e0   ;   rotn (:,:,:) = 0.e0 
    99          hdivb(:,:,:) = 0.e0   ;   hdivn(:,:,:) = 0.e0 
     96         !   before fields      !       now fields      
     97         sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp 
     98         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp 
     99         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp   
     100         rotb (:,:,:) = 0._wp   ;   rotn (:,:,:) = 0._wp 
     101         hdivb(:,:,:) = 0._wp   ;   hdivn(:,:,:) = 0._wp 
     102         ! 
     103         !                                       ! define e3u_b, e3v_b from e3t_b initialized in domzgr 
     104         CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 
    100105         ! 
    101106         IF( cp_cfg == 'eel' ) THEN 

  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r2779 r3008  
    9292      !!               un,vn   now horizontal velocity of next time-step 
    9393      !!---------------------------------------------------------------------- 
    94       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    9594      USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace 
    96       USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3 
    9795      ! 
    9896      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    9997      ! 
    10098      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     99      INTEGER  ::   iku, ikv     ! local integers 
    101100#if ! defined key_dynspg_flt 
    102101      REAL(wp) ::   z2dt         ! temporary scalar 
    103102#endif 
    104       REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars 
    105       REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      - 
    106       REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1 
     103      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars 
     104      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      - 
    107105      !!---------------------------------------------------------------------- 
    108  
    109       IF( wrk_in_use(2, 1,2,3) ) THEN 
    110          CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN 
    111       ENDIF 
    112106 
    113107      IF( kt == nit000 ) THEN 
    … …  
    238232         ELSE                             ! Variable volume ! 
    239233            !                             ! ================! 
    240             ! Before scale factor at t-points 
    241             ! ------------------------------- 
    242             DO jk = 1, jpkm1 
     234            ! 
     235            DO jk = 1, jpkm1                 ! Before scale factor at t-points 
    243236               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   & 
    244237                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     & 
    245                   &                         - 2.e0 * fse3t_n(:,:,jk)            ) 
    246             ENDDO 
    247             ! Add volume filter correction only at the first level of t-point scale factors 
    248             zec = atfp * rdt / rau0 
     238                  &                         - 2._wp * fse3t_n(:,:,jk)            ) 
     239            END DO 
     240            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors 
    249241            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    250             ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 
    251             zs_t  (:,:) =       e1t(:,:) * e2t(:,:) 
    252             zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) ) 
    253             zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) ) 
    254242            ! 
    255             IF( ln_dynadv_vec ) THEN 
    256                ! Before scale factor at (u/v)-points 
    257                ! ----------------------------------- 
    258                ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    259                DO jk = 1, jpkm1 
    260                   DO jj = 1, jpjm1 
    261                      DO ji = 1, jpim1 
    262                         zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    263                         zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    264                         zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    265                         fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    266                         fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
    267                      END DO 
    268                   END DO 
    269                END DO 
    270                ! lateral boundary conditions 
    271                CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) 
    272                CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
    273                ! Add initial scale factor to scale factor anomaly 
    274                fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
    275                fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 
    276                ! Leap-Frog - Asselin filter and swap: applied on velocity 
    277                ! ----------------------------------- 
    278                DO jk = 1, jpkm1 
    279                   DO jj = 1, jpj 
     243            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation) 
     244               ! 
     245               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b) 
     246               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 
     247               ! 
     248               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity 
     249                  DO jj = 1, jpj                      !                                                 -------- 
    280250                     DO ji = 1, jpi 
    281251                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    … …  
    290260               END DO 
    291261               ! 
    292             ELSE 
    293                ! Temporary filered scale factor at (u/v)-points (will become before scale factor) 
    294                !----------------------------------------------- 
    295                ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    296                DO jk = 1, jpkm1 
    297                   DO jj = 1, jpjm1 
    298                      DO ji = 1, jpim1 
    299                         zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    300                         zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    301                         zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    302                         ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    303                         ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
    304                      END DO 
    305                   END DO 
    306                END DO 
    307                ! lateral boundary conditions 
    308                CALL lbc_lnk( ze3u_f, 'U', 1. ) 
    309                CALL lbc_lnk( ze3v_f, 'V', 1. ) 
    310                ! Add initial scale factor to scale factor anomaly 
    311                ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 
    312                ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 
    313                ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 
    314                ! -----------------------------------             =========================== 
    315                DO jk = 1, jpkm1 
    316                   DO jj = 1, jpj 
    317                      DO ji = 1, jpim1 
     262            ELSE                             ! flux form (thickness weighted calulation) 
     263               ! 
     264               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b) 
     265               ! 
     266               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap:  
     267                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity 
     268                     DO ji = 1, jpim1                 !                              --------------------------- 
    318269                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 
    319270                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk) 
    … …  
    323274                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk) 
    324275                        ! 
    325                         zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
    326                         zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
     276                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
     277                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
    327278                        ! 
    328                         ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     279                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity 
    329280                        vb(ji,jj,jk) = zvf 
    330                         un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
     281                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua 
    331282                        vn(ji,jj,jk) = va(ji,jj,jk) 
    332283                     END DO 
    333284                  END DO 
    334285               END DO 
    335                fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor 
    336                fse3v_b(:,:,:) = ze3v_f(:,:,:) 
    337                CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions 
     286               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor 
     287               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
     288               CALL lbc_lnk( ub, 'U', -1. )                    ! lateral boundary conditions 
    338289               CALL lbc_lnk( vb, 'V', -1. ) 
    339290            ENDIF 
    … …  
    346297         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    347298      !  
    348       IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays') 
    349       ! 
    350299   END SUBROUTINE dyn_nxt 
    351300 

  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r2724 r3008  
    119119      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    120120      INTEGER  ::   icycle           ! local scalar 
    121       REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! local scalars 
     121      REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b     ! local scalars 
    122122      REAL(wp) ::   z1_8, zx1, zy1                   !   -      - 
    123123      REAL(wp) ::   z1_4, zx2, zy2                   !   -      - 
    … …  
    161161 
    162162      !                                   !* Local constant initialization 
    163       z2dt_b = 2.0 * rdt                                    ! baroclinic time step 
     163      z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                               ! baroclinic time step 
     164      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! baroclinic time step (starting with euler timestep) 
    164165      z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates 
    165166      z1_4 = 0.5 * 0.5 
    … …  
    195196               ! 
    196197#if defined key_vvl 
    197                ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk)  
    198                vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)    
     198!              ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk)  
     199!              vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)    
     200               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk)  
     201               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk) 
    199202#else 
    200203               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
    … …  
    278281      !                                             ! Remove barotropic contribution of bottom friction  
    279282      !                                             ! from the barotropic transport trend 
    280       zcoef = -1. / z2dt_b 
     283      zcoef = -1._wp * z1_2dt_b 
    281284# if defined key_vectopt_loop 
    282285      DO jj = 1, 1 
    … …  
    565568      !                                   !* update the general momentum trend 
    566569      DO jk=1,jpkm1 
    567          ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b 
    568          va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b 
     570         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 
     571         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 
    569572      END DO 
    570573      un_b  (:,:) =  zu_sum(:,:)  
Note: See TracChangeset for help on using the changeset viewer.