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 5038 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 – NEMO

Ignore:
Timestamp:
2015-01-20T15:26:13+01:00 (9 years ago)
Author:
jamesharle
Message:

Merging branch with HEAD of the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r4370 r5038  
    1818   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 
    1919   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
     20   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends 
    2021   !!------------------------------------------------------------------------- 
    2122   
     
    3435   USE bdydyn          ! ocean open boundary conditions 
    3536   USE bdyvol          ! ocean open boundary condition (bdy_vol routines) 
     37   USE trd_oce         ! trends: ocean variables 
     38   USE trddyn          ! trend manager: dynamics 
     39   USE trdken          ! trend manager: kinetic energy 
     40   ! 
    3641   USE in_out_manager  ! I/O manager 
     42   USE iom             ! I/O manager library 
    3743   USE lbclnk          ! lateral boundary condition (or mpp link) 
    3844   USE lib_mpp         ! MPP library 
    3945   USE wrk_nemo        ! Memory Allocation 
    4046   USE prtctl          ! Print control 
    41  
     47   USE timing          ! Timing 
    4248#if defined key_agrif 
    4349   USE agrif_opa_interp 
    4450#endif 
    45    USE timing          ! Timing 
    4651 
    4752   IMPLICIT NONE 
     
    7984      !!             at the local domain boundaries through lbc_lnk call, 
    8085      !!             at the one-way open boundaries (lk_bdy=T), 
    81       !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
     86      !!             at the AGRIF zoom   boundaries (lk_agrif=T) 
    8287      !! 
    8388      !!              * Apply the time filter applied and swap of the dynamics 
     
    99104      REAL(wp) ::   z2dt         ! temporary scalar 
    100105#endif 
    101       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars 
    102       REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      - 
    103       REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva 
    104       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f  
     106      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars 
     107      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
     108      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve 
     109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva  
    105110      !!---------------------------------------------------------------------- 
    106111      ! 
    107       IF( nn_timing == 1 )  CALL timing_start('dyn_nxt') 
    108       ! 
    109       CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
    110       IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva ) 
     112      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt') 
     113      ! 
     114      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva ) 
     115      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve ) 
    111116      ! 
    112117      IF( kt == nit000 ) THEN 
     
    152157 
    153158# if defined key_dynspg_ts 
     159!!gm IF ( lk_dynspg_ts ) THEN .... 
    154160      ! Ensure below that barotropic velocities match time splitting estimate 
    155161      ! Compute actual transport and replace it with ts estimate at "after" time step 
    156       zua(:,:) = 0._wp 
    157       zva(:,:) = 0._wp 
    158       DO jk = 1, jpkm1 
    159          zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    160          zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
     162      zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 
     163      zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 
     164      DO jk = 2, jpkm1 
     165         zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     166         zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
    161167      END DO 
    162168      DO jk = 1, jpkm1 
    163          ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
    164          va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
     169         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
     170         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
    165171      END DO 
    166172 
     
    175181         END DO   
    176182      ENDIF 
     183!!gm ENDIF 
    177184# endif 
    178185 
     
    195202# endif 
    196203#endif 
     204 
     205      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics 
     206         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step  
     207         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt 
     208         ! 
     209         !                                  ! Kinetic energy and Conversion 
     210         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt ) 
     211         ! 
     212         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends 
     213            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 
     214            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 
     215            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter 
     216            CALL iom_put( "vtrd_tot", zva ) 
     217         ENDIF 
     218         ! 
     219         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter 
     220         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the 
     221         !                                  !  computation of the asselin filter trends) 
     222      ENDIF 
    197223 
    198224      ! Time filter and swap of dynamics arrays 
     
    217243               DO jj = 1, jpj 
    218244                  DO ji = 1, jpi     
    219                      zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    220                      zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 
     245                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 
     246                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 
    221247                     ! 
    222248                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     
    301327            ! Revert "before" velocities to time split estimate 
    302328            ! Doing it here also means that asselin filter contribution is removed   
    303             zua(:,:) = 0._wp 
    304             zva(:,:) = 0._wp 
    305             DO jk = 1, jpkm1 
    306                zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
    307                zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
     329            zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 
     330            zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)     
     331            DO jk = 2, jpkm1 
     332               zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
     333               zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
    308334            END DO 
    309335            DO jk = 1, jpkm1 
    310                ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk) 
    311                vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk) 
     336               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk) 
     337               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk) 
    312338            END DO 
    313339         ENDIF 
     
    327353            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    328354         END DO 
    329          hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) ) 
    330          hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) ) 
     355         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 
     356         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 
    331357      ENDIF 
    332358      ! 
     
    335361      ! 
    336362      DO jk = 1, jpkm1 
    337 #if defined key_vectopt_loop 
    338          DO jj = 1, 1         !Vector opt. => forced unrolling 
    339             DO ji = 1, jpij 
    340 #else  
    341363         DO jj = 1, jpj 
    342364            DO ji = 1, jpi 
    343 #endif                   
    344365               un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 
    345366               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
     
    358379      ! 
    359380      ! 
     381 
     382      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum 
     383         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 
     384         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 
     385         CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 
     386      ENDIF 
     387      ! 
    360388      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    361389         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    362390      !  
    363       CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
    364       IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva ) 
     391      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva ) 
     392      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve ) 
    365393      ! 
    366394      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt') 
     
    370398   !!========================================================================= 
    371399END MODULE dynnxt 
    372  
Note: See TracChangeset for help on using the changeset viewer.