Changeset 2807


Ignore:
Timestamp:
2011-07-19T18:35:40+02:00 (9 years ago)
Author:
acc
Message:

Branch: dev_r2802_NOCS_vvlfix. Bugfix corrections to the calculation of fse3u_b and fse3v_b to address problems with vvl and partial steps. See comments added to ticket #812

Location:
branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r2779 r2807  
    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) ) 
     166      ! 
     167      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays') 
     168      ! 
     169   END SUBROUTINE dom_vvl 
     170 
     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      REAL(wp) ::   zvt          ! local scalars 
     191      !!---------------------------------------------------------------------- 
     192 
     193      IF( lwp .AND. kt == nit000 ) THEN 
     194         WRITE(numout,*) 
     195         WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization' 
     196         WRITE(numout,*) '~~~~~~~~~ ' 
     197         pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk) 
     198         pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk) 
     199      ENDIF 
     200       
     201      DO jk = 1, jpkm1           ! set the before scale factors at u- & v-points 
     202         DO jj = 2, jpjm1 
     203            DO ji = fs_2, fs_jpim1 
     204               zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj) 
     205               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) ) 
     206               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) ) 
    182207            END DO 
    183208         END DO 
    184209      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') 
    192       ! 
    193    END SUBROUTINE dom_vvl 
    194  
     210      IF( ln_zps ) THEN          ! minimum of the e3t at partial cell level 
     211         DO jj = 2, jpjm1 
     212            DO ji = fs_2, fs_jpim1 
     213               iku = mbku(ji,jj) 
     214               ikv = mbkv(ji,jj) 
     215               pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj  ,iku) )  
     216               pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji  ,jj+1,ikv) )  
     217            END DO 
     218         END DO 
     219      ENDIF 
     220      pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:)      ! anomaly to avoid zero along closed boundary/extra halos 
     221      pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:) 
     222      CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions 
     223      CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. ) 
     224      pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:)      ! recover the full scale factor 
     225      pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:) 
     226      ! 
     227   END SUBROUTINE dom_vvl_2 
     228    
    195229#else 
    196230   !!---------------------------------------------------------------------- 
     
    200234   SUBROUTINE dom_vvl 
    201235   END SUBROUTINE dom_vvl 
     236   SUBROUTINE dom_vvl_2 
     237   END SUBROUTINE dom_vvl_2 
    202238#endif 
    203239 
  • branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r2777 r2807  
    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_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r2779 r2807  
    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 
Note: See TracChangeset for help on using the changeset viewer.