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.
#1884 (Bug in diaar5.F90 (v3.6_STABLE & trunk)) – NEMO

Opened 7 years ago

Closed 7 years ago

Last modified 5 years ago

#1884 closed Bug (fixed)

Bug in diaar5.F90 (v3.6_STABLE & trunk)

Reported by: gm Owned by: timgraham
Priority: low Milestone:
Component: OCE Version: v3.6
Severity: Keywords:
Cc:

Description (last modified by nicolasmartin)

Context

key_diaar5 defined (diagnostics used for CMIP6)
When double diffusion parameterization is used (key_zdfddm), the diagnostic of the work done against stratification by vertical mixing ( iom_use( 'tnpeo' ) ) is wrong
In addition there is a out-of-bound : the jk-do loop should start at jk=2, not 1 ! In case of lk_zdfddm=F we can also start with jk=2 since act(:,:,1) is always zero.

This concerns the v3.6_STABLE and trunk versions.

Analysis

The portion of diaar5.F90 involved in the calculation is:

      IF( iom_use( 'tnpeo' )) THEN    
      ! Work done against stratification by vertical mixing
      ! Exclude points where rn2 is negative as convection kicks in here and
      ! work is not being done against stratification
          CALL wrk_alloc( jpi, jpj, zpe )
          zpe(:,:) = 0._wp
          IF( lk_zdfddm ) THEN
             DO ji=1,jpi
                DO jj=1,jpj
                   DO jk=1,jpk
                      zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
                         &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
                      !
                      zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
                      zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
                      !
                      zpe(ji, jj) = zpe(ji, jj) - MIN(0._wp, rn2(ji,jj,jk)) * &
                           &       grav * (  avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
                           &               - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )

                   ENDDO
                ENDDO
             ENDDO
          ELSE
             DO ji = 1, jpi
                DO jj = 1, jpj
                   DO jk = 1, jpk
                       zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
                   ENDDO
                ENDDO
             ENDDO
          ENDIF
          CALL lbc_lnk( zpe, 'T', 1._wp)         
          CALL iom_put( 'tnpeo', zpe )
          CALL wrk_dealloc( jpi, jpj, zpe )
      ENDIF

When lk_zdfddm = F, the calculation is OK and corresponds to a vertical sum of rau0 * avt * N2
When lk_zdfddm = T, the loop involved both a multiplication by rn2 (more precisely MIN( 0, rn2 ) ) AND by the vertical diffusive density flux taking into account avt and avs.
It is almost as if the vertical sum of rau0 * avt * N4 were computed (avt is not significantly different from avs in most of the world ocean...)

Fix

A solution is to replace the lk_zdfddm = True case by:

         IF( ln_zdfddm ) THEN
            DO jk = 2, jpk
               DO jj = 1, jpj
                  DO ji = 1, jpi
                     IF( rn2(ji,jj,jk) > 0._wp ) THEN
                        zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
                           &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
                        !
                        zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
                        zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
                        !
                        zpe(ji, jj) = zpe(ji, jj)            &
                           &        -  grav * (    avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
                           &                   - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
                     ENDIF
                  END DO
               END DO
             END DO
          ELSE
            DO jk = 1, jpk
               DO ji = 1, jpi
                  DO jj = 1, jpj
                     zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
                  END DO
               END DO
            END DO
         ENDIF

NB: I also invert the-do-loop order (i.e. from i-j-k to k-j-i) as it is faster due to the model 3D array structure

Commit History (6)

ChangesetAuthorTimeChangeLog
10870frrh2019-04-15T13:00:50+02:00

Apply correction to tnpeo calculation including reordering of k loop as per #1884. This version of the code reflects the vn 3.6 STABLE branch more than the vn 4.0 trunk. Although functionality is the same, here we include the lbc_lnk of the STABLE branch (not used in the trunk) and the workspace allocation of zpe as a pointer rather than the trunk version which defines zpe (and other workspace arrays) simply as normal allocatable arrays. i.e. I've opted to reflect the STABLE branch exactly rather than to transplant more optimal code from the trunk out of context.

8220lovato2017-06-26T18:17:34+02:00

3.6 stable: reinstate substitute variables in diaar5.F90 to compute tnpeo (#1884)

The substitute variables were replaced with the static ones in r8077, but it produces a compilation error when key_vvl is not used.

8084timgraham2017-05-26T14:04:22+02:00

Correction to commit for ticket #1884

8083timgraham2017-05-26T14:00:30+02:00

Correction to commit r8077 for ticket #1884

8078timgraham2017-05-26T11:16:21+02:00

#1884 fix in diaar5.F90 for calculation of tnpeo

8077timgraham2017-05-26T11:13:07+02:00

#1884 fix bug in diaar5 for calculation of tnpeo

Change History (6)

comment:1 Changed 7 years ago by gm

  • Description modified (diff)

comment:2 Changed 7 years ago by gm

  • Description modified (diff)

comment:3 Changed 7 years ago by nicolasmartin

  • Description modified (diff)

comment:4 Changed 7 years ago by clevy

  • Owner changed from nemo to timgraham

comment:5 Changed 7 years ago by timgraham

  • Resolution set to fixed
  • Status changed from new to closed

Fixed at revision r8077 of nemo_v3_6_STABLE and r8078 of the trunk

Last edited 7 years ago by timgraham (previous) (diff)

comment:6 Changed 5 years ago by frrh

In 10870:

Error: Failed to load processor CommitTicketReference
No macro or processor named 'CommitTicketReference' found
Note: See TracTickets for help on using tickets.