Skip to content
Snippets Groups Projects
Commit e46bd51c authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

correct a small bug in z restart

parent a0eab701
No related branches found
No related tags found
No related merge requests found
...@@ -133,8 +133,8 @@ CONTAINS ...@@ -133,8 +133,8 @@ CONTAINS
zjp1_cp = z_cp_stretched(z_idx_mapping(izg)+1) zjp1_cp = z_cp_stretched(z_idx_mapping(izg)+1)
zerotoone = (zi - zj_cp)/(zjp1_cp-zj_cp) ! weight for interpolation zerotoone = (zi - zj_cp)/(zjp1_cp-zj_cp) ! weight for interpolation
moments(ia,ipi,iji,iky,ikx,izi,:) = & moments(ia,ipi,iji,iky,ikx,izi,:) = &
zerotoone *moments_cp(iacp,ipcp,ijcp,iycp,ixcp,z_idx_mapping(izg))& (1._xp-zerotoone) *moments_cp(iacp,ipcp,ijcp,iycp,ixcp,z_idx_mapping(izg))&
+(1._xp-zerotoone)*moments_cp(iacp,ipcp,ijcp,iycp,ixcp,z_idx_mapping(izg+1)) + zerotoone *moments_cp(iacp,ipcp,ijcp,iycp,ixcp,z_idx_mapping(izg+1))
ELSE ! copy on all planes ELSE ! copy on all planes
moments(ia,ipi,iji,iky,ikx,izi,:) = moments_cp(iacp,ipcp,ijcp,iycp,ixcp,1) moments(ia,ipi,iji,iky,ikx,izi,:) = moments_cp(iacp,ipcp,ijcp,iycp,ixcp,1)
ENDIF ENDIF
...@@ -177,35 +177,45 @@ CONTAINS ...@@ -177,35 +177,45 @@ CONTAINS
REAL(xp), DIMENSION(Nz_cp+1) , INTENT(OUT) :: z_cp_stretched REAL(xp), DIMENSION(Nz_cp+1) , INTENT(OUT) :: z_cp_stretched
! local variables ! local variables
INTEGER :: iz_new, jz_cp INTEGER :: iz_new, jz_cp
REAL(xp):: zi, zj_cp, dz_s REAL(xp):: zi, zj_cp, dz_s, dz_new
LOGICAL :: in_interval LOGICAL :: in_interval
! stretched checkpoint grid interval ! stretched checkpoint grid interval
dz_s = (zarray_new(Nz_new)-zarray_new(1))/(Nz_cp-1) dz_s = (zarray_new(Nz_new)-zarray_new(1))/(Nz_cp-1)
! We loop over each new z grid points dz_new= (zarray_new(Nz_new)-zarray_new(1))/(Nz_new-1)
DO iz_new = 1,Nz_new IF(ABS(dz_s-dz_new) .LT. EPSILON(dz_s)) THEN
zi = zarray_new(iz_new) ! current position DO iz_new = 1,Nz_new
! Loop over the stretched checkpoint grid to find the right interval z_cp_stretched(iz_new) = zarray_new(iz_new)
zj_cp = zarray_new(1) ! init cp grid position z_idx_mapping(iz_new) = iz_new
jz_cp = 1 ! init cp grid index ENDDO
in_interval = .FALSE. ! flag to check if we stand in the interval z_cp_stretched(Nz_new+1) = zarray_new(1)
DO WHILE (.NOT. in_interval) z_idx_mapping(Nz_new+1) = 1
in_interval = (zi .GE. zj_cp) .AND. (zi .LT. zj_cp+dz_s) ELSE ! Strech the grid
! Increment ! We loop over each new z grid points
zj_cp = zj_cp + dz_s DO iz_new = 1,Nz_new
jz_cp = jz_cp + 1 zi = zarray_new(iz_new) ! current position
IF(jz_cp .GT. Nz_cp+1) ERROR STOP "STOP: could not adapt grid .." ! Loop over the stretched checkpoint grid to find the right interval
ENDDO ! per construction the while loop should always top zj_cp = zarray_new(1) ! init cp grid position
z_idx_mapping(iz_new) = jz_cp-1 ! The last index was one too much so we store the one before jz_cp = 1 ! init cp grid index
ENDDO in_interval = .FALSE. ! flag to check if we stand in the interval
! we build explicitly the stretched cp grid for output and double check DO WHILE (.NOT. in_interval)
DO jz_cp = 1,Nz_cp in_interval = (zi .GE. zj_cp) .AND. (zi .LT. zj_cp+dz_s)
z_cp_stretched(jz_cp) = zarray_new(1) + (jz_cp-1)*dz_s ! Increment
ENDDO zj_cp = zj_cp + dz_s
! Periodicity jz_cp = jz_cp + 1
z_cp_stretched(Nz_cp+1) = z_cp_stretched(1) IF(jz_cp .GT. Nz_cp+1) ERROR STOP "STOP: could not adapt grid .."
z_idx_mapping (Nz_new+1) = z_idx_mapping (1) ENDDO ! per construction the while loop should always top
! Check that the start and the end of the grids are the same z_idx_mapping(iz_new) = jz_cp-1 ! The last index was one too much so we store the one before
IF(.NOT.(abs(z_cp_stretched(1)-zarray_new(1)) .LT. 1e-3).AND.(abs(z_cp_stretched(Nz_cp)-zarray_new(Nz_new)).LT.1e-3)) & ENDDO
ERROR STOP "Failed to stretch the cp grid" ! we build explicitly the stretched cp grid for output and double check
DO jz_cp = 1,Nz_cp
z_cp_stretched(jz_cp) = zarray_new(1) + (jz_cp-1)*dz_s
ENDDO
! Periodicity
z_cp_stretched(Nz_cp+1) = z_cp_stretched(1)
z_idx_mapping (Nz_new+1) = z_idx_mapping(1)
! Check that the start and the end of the grids are the same
IF(.NOT.(abs(z_cp_stretched(1)-zarray_new(1)) .LT. 1e-3).AND.(abs(z_cp_stretched(Nz_cp)-zarray_new(Nz_new)).LT.1e-3)) &
ERROR STOP "Failed to stretch the cp grid"
ENDIF
END SUBROUTINE z_grid_mapping END SUBROUTINE z_grid_mapping
END MODULE restarts END MODULE restarts
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment