Opened 3 years ago

Last modified 7 weeks ago

#1584 assigned Bug

utrd_zdf and vtrd_zdf are velocities minus momentum trends

Reported by:chrenkl Component: OPA
Version: release-3.6 Severity: minor
Keywords: TRD diagnostics Cc: Gurvan Madec
Management
Owned by:Gurvan Madec Milestone:
Priority: high

Description (last modified by Nicolas MARTIN)

In the subroutine for vertical diffusion using a implicit time-stepping (dyn_zdf_imp in dynzdf_imp.F90), the momentum trends ua and va are integrated in time and become the after velocities (lines 122-136). Thus, in lines 90-95 of dynzdf.F90, when the momentum trend diagnostics are supposed to be computed, momentum trends (ztrdu, ztrdv) are subtracted from velocities (ua, va).

Therefore,

90   IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
91      ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
92      ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
93      CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
94      CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
95   ENDIF

should be changed to

90   IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
91      ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
92      ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdu(:,:,:)
93      CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
94      CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
95   ENDIF

as it is done in dynspg.F90. I have not checked whether this is also the case for vertical diffusion with explicit time-stepping.

Commit History (2)

ChangesetAuthorTimeChangeLog
6132jchanut2015-12-18T14:54:06+01:00

Correct momentum vertical diffusion tendencies, #1584

6131jchanut2015-12-18T14:52:13+01:00

Correct momentum vertical diffusion tendencies, #1584

Attachments (1)

dynzdf.F90 (8.8 KB) - added by chrenkl 5 months ago.
Modified dynzdf.F90

Download all attachments as: .zip

Change History (15)

comment:1 Changed 20 months ago by Nicolas MARTIN

  • Keywords trends added; trend removed

comment:2 Changed 20 months ago by Nicolas MARTIN

  • Keywords nemo_v3_6* added

comment:3 Changed 20 months ago by Nicolas MARTIN

  • Keywords TRD nemo_v3_6_STABLE added; trends removed

comment:4 Changed 20 months ago by Nicolas MARTIN

  • Keywords nemo_v3_6_STABLE removed

comment:5 Changed 9 months ago by Claire Levy

  • Owner changed from nemo to lovato

comment:6 Changed 9 months ago by Tomas Lovato

  • Owner changed from lovato to gm

comment:7 Changed 5 months ago by Claire Levy

  • Cc gm added
  • Status changed from new to assigned

comment:8 Changed 5 months ago by Nicolas MARTIN

  • Priority changed from minor to major

comment:9 Changed 5 months ago by Nicolas MARTIN

  • Description modified (diff)

comment:10 Changed 5 months ago by Claire Levy

  • Priority changed from major to critical

comment:11 Changed 5 months ago by chrenkl

It seems like my initial solution did not quite work. I have since modified the code (see attachment). This might not be the most straightforward solution, but it seems to work.

Changed 5 months ago by chrenkl

Modified dynzdf.F90

comment:12 Changed 3 months ago by nemo

  • Keywords release-3.6* added; nemo_v3_6* removed

comment:13 Changed 3 months ago by nemo

  • Keywords release-3.6* removed

comment:14 Changed 7 weeks ago by Nicolas MARTIN

  • Severity set to minor

The status of this ticket is not clear.

2 commits (r6131 & r6132) have been linked to it dating back to the end of 2015 but the ticket had not been closed. Then the ZDF routines have been vastly re-factored in 2017 before merging with the trunk (see #1883).
Anyway the initial fix is already applied in the new routine:

Line/branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
384      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
385         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
386         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
387         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
388         DEALLOCATE( ztrdu, ztrdv ) 
389      ENDIF


But a decision has to be made regarding the last fix suggested by Christoph (extract of dynzdf.F90 attached)

   SUBROUTINE dyn_zdf( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_zdf  ***
      !!
      !! ** Purpose :   compute the vertical ocean dynamics physics.
      !!---------------------------------------------------------------------
      !!
      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
      !
      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
! CR 20161028 --- begin ---
      INTEGER                             ::  jk   ! dummy loop indices
      REAL(wp)                            ::  z1_2dt 
      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zua, zva 
! CR 20161028 --- end ---
      !!---------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf')
! CR 20161028 --- begin ---
!       CALL wrk_alloc( jpi,jpj,jpk, zua, zva )
!       CALL wrk_alloc( jpi,jpj, zue, zve )
! CR 20161028 --- end ---
      !
      !                                          ! set time step
      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdtra (restart with Euler time stepping)
      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdttra (leapfrog)
      ENDIF

      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
         ztrdu(:,:,:) = ua(:,:,:)
         ztrdv(:,:,:) = va(:,:,:)
      ENDIF

      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
      !
      CASE ( 0 )   ;   CALL dyn_zdf_exp( kt, r2dt )      ! explicit scheme
      CASE ( 1 )   ;   CALL dyn_zdf_imp( kt, r2dt )      ! implicit scheme
      !
      CASE ( -1 )                                        ! esopa: test all possibility with control print
                       CALL dyn_zdf_exp( kt, r2dt )
                       CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf0 - Ua: ', mask1=umask,               &
                          &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
                       CALL dyn_zdf_imp( kt, r2dt )
                       CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf1 - Ua: ', mask1=umask,               &
                          &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
      END SELECT

      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
! CR 20161028         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
! CR 20161028         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
! CR 20161028 --- begin ---
      ! In dyn_zdf_imp, ua & va become velocities and are not trends anymore. Thus, we have to back-calculate
      ! in order to obtain the trend due to vertical diffusion.
      ! Before, replace barotropic component: (copied from dynnxt.F90))
      ! Ensure below that barotropic velocities match time splitting estimate
      ! Compute actual transport and replace it with ts estimate at "after" time step
         CALL wrk_alloc( jpi,jpj,jpk, zua, zva )
         CALL wrk_alloc( jpi,jpj, zue, zve )
         zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
         zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
         DO jk = 2, jpkm1
            zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
            zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
         END DO
         DO jk = 1, jpkm1
            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
         END DO
         !
         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt     
         !
         zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt - ztrdu
         zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt - ztrdv
         CALL trd_dyn( zua, zva, jpdyn_zdf, kt )
         CALL wrk_dealloc( jpi,jpj,jpk, zua, zva )
         CALL wrk_dealloc( jpi,jpj, zue, zve )
! CR 20161028 --- end ---
! CR 20161031         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
      ENDIF
      !                                          ! print mean trends (used for debugging)
      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
            &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
      !
      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf')
      !
   END SUBROUTINE dyn_zdf
Note: See TracTickets for help on using tickets.