| Version 15 (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.
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.
r14549 Allow RK3 time-stepping with 2D mode.
Prepare forcings and barotropic 2D fields. dynspg_ts remains for 2D mode integration. dyn_vor_RK3 only computes 2D relative vorticity and metric term from 3D to 2D.
stp_2D is called by stprk3 in single first.
OCE
|--DYN
|-- dynspg_ts.F90
#remove k_only_ADV
PUBLIC dyn_drg_init ! called by rk3ssh
! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
IF( kt == nit000 ) THEN
IF( .NOT.ln_bt_fw .OR. ln_bt_av ) CALL ctl_stop( 'dyn_spg_ts: RK3 requires ln_bt_fw=T AND ln_bt_av=F')
ENDIF
! ! set values computed in RK3_ssh
zssh_frc(:,:) = sshe_rhs(:,:)
zu_frc(:,:) = Ue_rhs(:,:)
zv_frc(:,:) = Ve_rhs(:,:)
zCdU_u (:,:) = CdU_u (:,:)
zCdU_v (:,:) = CdU_v (:,:)
!
IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient
! Phase 3. update the general trend with the barotropic trend
IF(.NOT.ln_bt_av ) THEN !* Update Kaa barotropic external mode
uu_b(:,:,Kaa) = ua_e (:,:)
pvv_b(:,:,Kaa) = va_e (:,:)
pssh (:,:,Kaa) = ssha_e(:,:)
ENDIF
un_adv...
ubar...
|-- dynspg.F90
#remove k_only_ADV
|-- dynvor.F90
SUBROUTINE dyn_vor_3D( kt, Kmm, puu, pvv, Krhs )
|-- dynzdf.F90
zDt_2 = rDt * 0.5_wp # small cosmetic optim
|-- stprk3_stg.F90
|-- stp2d.F90
|-- stprk3.F90
179 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
180 ! RK3 : single first external mode computation
181 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
182
183 CALL stp_2D( kstp, Nbb, Nbb, Naa, Nrhs ) ! out: ssh, (uu_b,vv_b) and (un_adv,vn_adv) at Naa
184
185
186 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
187 ! RK3 time integration
188 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
189
190 ! Stage 1 :
191 CALL stp_RK3_stg( 1, kstp, Nbb, Nbb, Nrhs, Naa )
192 !
193 Nrhs = Nnn ; Nnn = Naa ; Naa = Nrhs ! Swap: Nbb unchanged, Nnn <==> Naa
194 !
195 ! Stage 2 :
196 CALL stp_RK3_stg( 2, kstp, Nbb, Nnn, Nrhs, Naa )
197 !
198 Nrhs = Nnn ; Nnn = Naa ; Naa = Nrhs ! Swap: Nbb unchanged, Nnn <==> Naa
199 !
200 ! Stage 3 :
201 CALL stp_RK3_stg( 3, kstp, Nbb, Nnn, Nrhs, Naa )
202 !
203 Nrhs = Nbb ; Nbb = Naa ; Naa = Nrhs ! Swap: Nnn unchanged, Nbb <==> Naa
Documentation updates
...
Preview
...
Tests
...
Review
...
