| 1 | subroutine time_step |
|---|
| 2 | #define DEBUG_LEVEL 1 |
|---|
| 3 | ! Routine to take an explicit time step via Euler, |
|---|
| 4 | ! 2nd or third order Adams Bashforth method. The differential |
|---|
| 5 | ! equations are as follows: |
|---|
| 6 | ! |
|---|
| 7 | ! d/dt u = rhs1 |
|---|
| 8 | ! d/dt v = rhs2 |
|---|
| 9 | ! d/dt w = rhs3 |
|---|
| 10 | ! d/dt s1 = rhs4 |
|---|
| 11 | ! d/dt s2 = rhs5 |
|---|
| 12 | ! |
|---|
| 13 | ! notation: |
|---|
| 14 | ! MM0 integer index identifying the time t storage locations |
|---|
| 15 | ! MM1 integer index identifying the time t-dt storage locations |
|---|
| 16 | ! MM2 integer index identifying the time t-2*dt storage locations |
|---|
| 17 | ! rhs1(:,:,:,MM0),rhs1(:,:,:,MM1),rhs1(:,:,:,MM2) rhs1 at t,t-dt,t-2*dt |
|---|
| 18 | ! etc for rhs2,rhs3,rhs4,rhs5 |
|---|
| 19 | ! dt time increment |
|---|
| 20 | ! |
|---|
| 21 | ! outputs: |
|---|
| 22 | ! u_star(1:nx,1:ny,1:nz) calculated values of u at time t+dt |
|---|
| 23 | ! v_star(1:nx,1:ny,1:nz) calculated values of v at time t+dt |
|---|
| 24 | ! w_star(1:nx,1:ny,1:nz) calculated values of w at time t+dt |
|---|
| 25 | ! s1(1:nx,1:ny,1:nz) calculated values of s1 at time t+dt (overwrite) |
|---|
| 26 | ! s2(1:nx,1:ny,1:nz) calculated values of s2 at time t+dt (overwrite) |
|---|
| 27 | |
|---|
| 28 | use grid_info, only: dt |
|---|
| 29 | use dependent_variables |
|---|
| 30 | use intermediate_variables |
|---|
| 31 | use rhs_variables |
|---|
| 32 | use counters_flags_etc |
|---|
| 33 | use mpi_parameters |
|---|
| 34 | implicit none |
|---|
| 35 | |
|---|
| 36 | include '../input/problem_size.h' |
|---|
| 37 | real dt2,dt12,dt24 |
|---|
| 38 | #if DEBUG_LEVEL >= 1 |
|---|
| 39 | if(myid==0) write(0,*) 'hello world from time_step' |
|---|
| 40 | #endif |
|---|
| 41 | |
|---|
| 42 | dt2=dt/2. |
|---|
| 43 | dt12=dt/12. |
|---|
| 44 | dt24=dt/24. |
|---|
| 45 | |
|---|
| 46 | if( step_flag == 'EULER') then |
|---|
| 47 | u_star(:,:,:)= u(:,:,:)+dt*rhs1(:,:,:,MM0) |
|---|
| 48 | v_star(:,:,:)= v(:,:,:)+dt*rhs2(:,:,:,MM0) |
|---|
| 49 | w_star(:,:,:)= w(:,:,:)+dt*rhs3(:,:,:,MM0) |
|---|
| 50 | s1(:,:,:) =s1(:,:,:)+dt*rhs4(:,:,:,MM0) |
|---|
| 51 | if(num_scalars>1) s2(:,:,:)=s2(:,:,:)+dt*rhs5(:,:,:,MM0) |
|---|
| 52 | step_flag='AB2' |
|---|
| 53 | elseif( step_flag == 'AB2' ) then |
|---|
| 54 | u_star(:,:,:)= u(:,:,:)+dt2*(3.*rhs1(:,:,:,MM0)-rhs1(:,:,:,MM1)) |
|---|
| 55 | v_star(:,:,:)= v(:,:,:)+dt2*(3.*rhs2(:,:,:,MM0)-rhs2(:,:,:,MM1)) |
|---|
| 56 | w_star(:,:,:)= w(:,:,:)+dt2*(3.*rhs3(:,:,:,MM0)-rhs3(:,:,:,MM1)) |
|---|
| 57 | s1(:,:,:) =s1(:,:,:)+dt2*(3.*rhs4(:,:,:,MM0)-rhs4(:,:,:,MM1)) |
|---|
| 58 | if(num_scalars>1) s2(:,:,:)=s2(:,:,:)+dt2*(3.*rhs5(:,:,:,MM0)-rhs5(:,:,:,MM1)) |
|---|
| 59 | if( AB_order >= 3 ) step_flag='AB3' |
|---|
| 60 | elseif( step_flag == 'AB3' ) then |
|---|
| 61 | u_star= u+dt12*(23.*rhs1(:,:,:,MM0)-16.*rhs1(:,:,:,MM1)+5.*rhs1(:,:,:,MM2)) |
|---|
| 62 | v_star= v+dt12*(23.*rhs2(:,:,:,MM0)-16.*rhs2(:,:,:,MM1)+5.*rhs2(:,:,:,MM2)) |
|---|
| 63 | w_star= w+dt12*(23.*rhs3(:,:,:,MM0)-16.*rhs3(:,:,:,MM1)+5.*rhs3(:,:,:,MM2)) |
|---|
| 64 | s1 =s1+dt12*(23.*rhs4(:,:,:,MM0)-16.*rhs4(:,:,:,MM1)+5.*rhs4(:,:,:,MM2)) |
|---|
| 65 | if(num_scalars>1) s2=s2+dt12*(23.*rhs5(:,:,:,MM0)-16.*rhs5(:,:,:,MM1)+5.*rhs5(:,:,:,MM2)) |
|---|
| 66 | if( AB_order == 4 ) step_flag='AB4' |
|---|
| 67 | elseif( step_flag == 'AB4' ) then |
|---|
| 68 | u_star= u+dt24*(55.*rhs1(:,:,:,MM0)-59.*rhs1(:,:,:,MM1)+37.*rhs1(:,:,:,MM2)-9.*rhs1(:,:,:,MM3)) |
|---|
| 69 | v_star= v+dt24*(55.*rhs2(:,:,:,MM0)-59.*rhs2(:,:,:,MM1)+37.*rhs2(:,:,:,MM2)-9.*rhs2(:,:,:,MM3)) |
|---|
| 70 | w_star= w+dt24*(55.*rhs3(:,:,:,MM0)-59.*rhs3(:,:,:,MM1)+37.*rhs3(:,:,:,MM2)-9.*rhs3(:,:,:,MM3)) |
|---|
| 71 | s1 =s1+dt24*(55.*rhs4(:,:,:,MM0)-59.*rhs4(:,:,:,MM1)+37.*rhs4(:,:,:,MM2)-9.*rhs4(:,:,:,MM3)) |
|---|
| 72 | if(num_scalars>1) s2=s2+dt24*(55.*rhs5(:,:,:,MM0)-59.*rhs5(:,:,:,MM1)+37.*rhs5(:,:,:,MM2)-9.*rhs5(:,:,:,MM3)) |
|---|
| 73 | endif |
|---|
| 74 | |
|---|
| 75 | end subroutine time_step |
|---|
