Name and subject of the action

Last edition: Wikinfo(changed_ts)? by Wikinfo(changed_by)?

The PI is responsible to closely follow the progress of the action, and especially to contact NEMO project manager if the delay on preview (or review) are longer than the 2 weeks expected.

1. Summary
2. Preview
3. Tests
4. Review

Summary

Action	RK3 stage 1
PI(S)	Gurvan et Sibylle
Digest	Run a GYRE configuration with new RK3 scheme
Dependencies	If any
Branch	source:/NEMO/branches/2021/dev_r14318_RK3_stage1
Previewer(s)	Gurvan
Reviewer(s)	Names
Ticket	#2605

Description

Error: Failed to load processor box

No macro or processor named 'box' found

RK3 time stepping implementation for NEMO includes at this stage dynamic and active tracers implementation, time spitting single first with 2D mode integration.

...

Implementation

Error: Failed to load processor box

No macro or processor named 'box' found

RK3 implementation is splitted up into :

• code preparation
• dynamic and active tracers (barocline)
• vertical physics (TKE) ?
• barotropic mode (barotrope)
• mass forcing
• passive tracers

Code preparation In order to preserve constancy property velocity for momentum and active tracers must be the same. Advection routines in flux form are modified to take (u,v,w) as an input argument. In order to use advection routines for the barotropic mode we need the possibility to de-activate vertical advection computation. Advection routines in flux and vector form are modified to take an optional argument (no_zad) to do so.

Barocline part For sake of simplicity we started to implement RK3 regarding a GYRE configuration validation with no barotrope mode (ssh, uu_b, un_adv are set to zero at each time step). Forcing have been removed except winds and heat flux. key_qco is active and vertical physics is modeled as constant with high viscosity coefficients.

• Prepare routines
  • Change eos divhor and sshwzv interface.
• Add RK3 time stepping routines
  • rk3stg deals with time integration at N+1/3, N+1/2 and N+1
  • stprk3 orchestrates

Barotrope part In order to validate 2D mode implementation we remove above zero forcing for barotropic variables mass forcing remains to zero.

• Prepare routines
  • Change dynadv, dynvor, dynspg_ts
• Add RK3 2D mode time stepping routines
  • rk3ssh prepare 2D forcing, get dynamics 2D RHS from 3D trends, integrate 2D mode

...

Implementation details : Code preparation

r14418 Allow an advective velocity to be passed as an argument.
3D velocity can be a pointer.

OCE
 |-- oce.F90
     REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:), TARGET   ::   uu   ,  vv     !: horizontal velocities        [m/s] 
     REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  , TARGET   ::   ww             !: vertical velocity            [m/s]

3D velocity added as an input argument of advective routines passed through dyn_adv

OCE
 |--DYN
     |-- dynadv.F90
         SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
         ...
         CALL dyn_adv_cen2( kt                , Kmm, puu, pvv, Krhs, pau, pav, paw )   ! 2nd order centered scheme 
         CALL dyn_adv_ubs ( kt           , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )   ! 3rd order UBS      scheme (UP3) 
     |-- dynadv_cen2.F90
         SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw )  
         ...       
         IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww) 
            zptu => pau(:,:,:) 
            ...
         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk) 
     |-- dynadv_ubs.F90
         SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw ) 
         ...
         IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww) 
            zptu => pau(:,:,:) 
            ...
         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk) 
 |--TRA  
     |-- traadv.F90
         SUBROUTINE tra_adv( kt, Kbb, Kmm, pts, Krhs, pau, pav, paw ) 
         ...
         IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww) 
            zptu => pau(:,:,:) 
            ...
         zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * ( zptu(ji,jj,jk) + usd(ji,jj,jk) ) 

Finally this new structure is used in step and tested with usual velocities

OCE
 |-- stpmlf.F90
     REAL(wp), TARGET     , DIMENSION(jpi,jpj,jpk) ::   zau, zav, zaw   ! advective velocity 
     ...
     zau(:,:,:) = uu(:,:,:,Nnn)        !!st   patch for MLF will be computed in RK3
     ...
     CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs, zau, zav, zaw )  ! advection (VF or FF)   ==> RHS
     ...
     CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs, zau, zav, zaw )  ! hor. + vert. advection ==> RHS 

Results should be exactly the same as the ones from from the trunk. It was not the case for an OVERFLOW. The use of ln_wAimp=T changes ww at the truncature in diawri.F90, and that produces a small error. This has been corrected.

r14428 Allow vertical advection to be de-activated with an optionnal input argument : no_zad.

3D velocity added as an input argument of advective routines passed through dyn_adv

OCE
 |--DYN
     |-- dynadv.F90
         SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
         ...
         CALL dyn_adv_cen2( kt                , Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )   ! 2nd order centered scheme 
         CALL dyn_adv_ubs ( kt           , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )   ! 3rd order UBS      scheme (UP3) 
     |-- dynadv_cen2.F90
         SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )  
         ...       
         IF( PRESENT( no_zad ) ) THEN  !==  No vertical advection  ==!   (except if linear free surface) 
            IF( ln_linssh ) THEN                ! linear free surface: advection through the surface z=0 
               DO_2D( 0, 0, 0, 0 ) 
                  zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm) 
                  zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm) 
                  puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj)   & 
                     &                                        / e3u(ji,jj,1,Kmm) 
                  pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj)   & 
                     &                                        / e3v(ji,jj,1,Kmm) 
               END_2D 
            ENDIF 
         ! 
         ELSE                          !==  Vertical advection  ==! 
            ...
     |-- dynadv_ubs.F90
         SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) 
         ...
         IF( PRESENT( no_zad ) ) THEN  !==  No vertical advection  ==!   (except if linear free surface) 
            IF( ln_linssh ) THEN                ! linear free surface: advection through the surface z=0 
               DO_2D( 0, 0, 0, 0 ) 
                  zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
                  zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm) 
                  puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj)   & 
                     &                                        / e3u(ji,jj,1,Kmm) 
                  pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj)   & 
                     &                                        / e3v(ji,jj,1,Kmm) 
               END_2D 
            ENDIF 
         ! 
         ELSE                          !==  Vertical advection  ==! 

Gurvan added a loop optimisation for dynzad.F90

OCE
 |--DYN
     |-- dynzad.F90
All the loops are now gather in a single one.  

Implementation details : barocline processing

r14547 Allow RK3 time-stepping with 2D mode damped.
div_hor interface and sshwzv interface have been changed accordingly for RK3. eos also changed in order to avoid gdep to be used as an input argument in the key_qco framework.

OCE
 |--DYN
     |-- divhor.F90
         SUBROUTINE div_hor_RK3( kt, Kbb, Kmm, puu, pvv, pe3divUh ) 
     |-- sshwzv.F90
         SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww )
 |--TRA
     |-- eosbn2.F90
         SUBROUTINE eos_insitu_New( pts, Knn, prd ) 

Time step no longer need to be doubled. rk3 routines are added to the code and stprk3 is called through nemogcm when key_RK3 is active.

OCE
 |--DOM
     |-- domain.F90
         #if defined key_RK3 
              rDt   =         rn_Dt 
 	      r1_Dt = 1._wp / rDt
          ...
 |-- nemogcm.F90
     # if defined key_RK3 
             USE stprk3 
     ...
 |-- stprk3.F90
 |-- stprk3-stg.F90

Has been tested and validated against an modified leap frog GYRE in the same configuration with the same namelist.

Documentation updates

Error: Failed to load processor box

No macro or processor named 'box' found

...

Preview

Error: Failed to load processor box

No macro or processor named 'box' found

...

Tests

Error: Failed to load processor box

No macro or processor named 'box' found

...

Review

Error: Failed to load processor box

No macro or processor named 'box' found

...