Version 11 (modified by techene, 3 years ago) (diff) |
---|
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.
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
RK3 time stepping implementation for NEMO includes at this stage dynamic and active tracers implementation, time spitting single first with 2D mode integration.
...
Implementation
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.
Documentation updates
...
Preview
...
Tests
...
Review
...