ORCA2 LIM3 (new options) not restartable

[clevy]

After adding key_vvl and key_dynspg_ts, the ORCA2_LIM3 configuration tested in SETTE appears not to be restartable any more.
If it helps, in SETTE test the first time step after restart (=76) is the same, and changes appear afterwards (from time step 77 onwards)
Seems to appear when adding the vvl option.

[timgraham]

I get the same behaviour on our IBM. It looks like the problem might arise in SBC. I'm looking into it at the moment.

[clem]

The problem might come from the first ocean level fse3t which is used in limthd.F90 and which must be read from the restart file. It looks good to me in case ice is not embedded. I am not sure otherwise.

[timgraham]

I've managed to get ORCA2LIM3 to be restartable by resetting uslp and vslp to zero at about line 235 of ldfslp.F90 as follows:

         !
         !                                            !* horizontal Shapiro filter
         uslp(:,:,:) = 0._wp
         vslp(:,:,:) = 0._wp
         DO jk = 2, jpkm1
            DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
               DO ji = 2, jpim1
                  uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
                     &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
                     &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
                     &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
                     &                       + 4.*  zwz(ji  ,jj  ,jk)                       )

This works and doesn't change results (in ORCA2LIM2 or ORCA2LIM3 at least) but at the moment I don't understand why it's necessary in the VVL case. Based on the code that follows it must be a difference in top level of uslp and vslp after a restart.

[timgraham]

The problem arises in ldf_slp_init within the section following this if statement:
IF( ln_traldf_hor .OR. ln_dynldf_hor ) THEN
at line 849 (in ldf_slp_init). This is included for the case of horizontal mixing in S-coordinates but is being triggered in this case because the ORCA2 reference cases use ln_dynldf_hor=True. It didn't cause a problem without vvl because fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) = 0

Possible fixes are
1) to reset uslp and vslp to zero at each time step (as above) or to add a check for S coordinates to this
2) Change the if statement in ldf_slp_init to
IF( ln_sco .AND. ( ln_traldf_hor .OR. ln_dynldf_hor ) ) THEN
so that this bit of the code is only called in the S coordinate case.

However, should the Z coordinate VVL case be treated in a similar way to S coordinates where ln_dynldf_hor is true horizontal diffusion (and requires key_ldfslp) and ln_dynldf_level is diffusion on model levels?

[timgraham]

I've applied fix 2 (see above) at revision 5016.

[timgraham]

This didn't solve the problem as I had expected in a full length test so I'm reopening the ticket.

[timgraham]

The fix above is required but there is a 2nd problem somewhere in zdftmx. If I compile without key_zdftmx then ORCA2LIM3 is restartable.

[timgraham]

The problem is in zdf_tmx_init. The variable az_tmx is calculates in the following section of code using fsdepw. In the vvl case fsdepw is different at the start of the LONG run compared to when it is read in from the restart file at the start of the SHORT run leading to the restartability test failure.
In zdf_tmx_init:

      ! Vertical structure (az_tmx)
      DO jj = 1, jpj                ! part independent of the level
         DO ji = 1, jpi
            zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
            zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
         END DO
      END DO
      DO jk= 1, jpk                 ! complete with the level-dependent part
         DO jj = 1, jpj
            DO ji = 1, jpi
               az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-fsdepw(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
            END DO
         END DO
      END DO

This can be fixed by changing fsdepw to gdepw:

      ! Vertical structure (az_tmx)
      DO jj = 1, jpj                ! part independent of the level
         DO ji = 1, jpi
            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
            zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
         END DO
      END DO
      DO jk= 1, jpk                 ! complete with the level-dependent part
         DO jj = 1, jpj
            DO ji = 1, jpi
               az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-gdepw_0(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
            END DO
         END DO
      END DO

This change is applied at r5021.

As with the earlier fix there is a question about whether this fix is correct in the long term. Should az_tmx actually be updated at every timestep in the VVL case? I'll put this in a new ticket.