MODULE tools_brice !----------------------------------------------------------- ! ! to make it we use a 4th order polynomial interpolation ! ! Created by Brice Lemaire on 12/2009. ! !----------------------------------------------------------- USE agrif_types ! IMPLICIT NONE PUBLIC ! INTEGER :: nxGmix, nyGmix INTEGER :: nxG1, nyG1, nxG2, nyG2, nxG3, nyG3, nxG4, nyG4 INTEGER :: nxG1mix, nyG1mix INTEGER :: nxG2mix, nyG2mix INTEGER :: nxG3mix, nyG3mix INTEGER :: nxG4mix, nyG4mix INTEGER :: nx_fine, ny_fine, nx_coarse, ny_coarse INTEGER :: i, j, i_min, j_min INTEGER :: k, m ! PUBLIC mixed_grid PRIVATE mixed_G1, mixed_G2, mixed_G4 ! INTERFACE mixed_grid MODULE PROCEDURE mixed_G1, mixed_G2, mixed_G4 END INTERFACE ! CONTAINS ! !******************************************************** ! SUBROUTINE mixed grid * ! * ! mix from four grids (U,V,F,T) to one grid * ! * ! CALLED from create_coordinates.f90 * !******************************************************** ! SUBROUTINE mixed_G1(G1,Gmix) ! TYPE(Coordinates), INTENT(IN) :: G1 TYPE(mixed_coordinates), INTENT(INOUT) :: Gmix ! WRITE(*,*) '' WRITE(*,*) '*** ROUTINE mixed_G1 ***' WRITE(*,*) '' ! nxG1 = (imax+1) - (imin-1) + 1 !(+1) rajout de 2 bandes fantmes nyG1 = (jmax+1) - (jmin-1) + 1 ! WRITE(*,*) '' WRITE(*,*) '*** CHECKING SIZE OF INPUT GRID***' WRITE(*,*) nxG1, 'x', nyG1 WRITE(*,*) '' ! nx_coarse = nxG1 ny_coarse = nyG1 ! !!!Calculate size of mixed grid (nx x ny) nxGmix = (nx_coarse) * 2 !nbre de pts connus (T,U,V,F) suivant x nxGmix = nxGmix + (rho-1)*(nxGmix-1) !nbre de points interpoler ! nyGmix = (ny_coarse) * 2 !nbre de pts connus (T,U,V,F) suivant y nyGmix = nyGmix + (rho-1)*(nyGmix-1) !nbre de points interpoler ! WRITE(*,*) '' WRITE(*,*) '*** SIZE OF MIXED GRID ***' WRITE(*,*) nxGmix, ' x ', nyGmix WRITE(*,*) '' ! nx_fine = (nx_coarse-2)*rho + 1 ny_fine = (ny_coarse-2)*rho + 1 ! WRITE(*,*) '' WRITE(*,*) '*** SIZE OF FINE GRID ***' WRITE(*,*) nx_fine, ' x ', ny_fine WRITE(*,*) '' ! !!!Allocate mixed grid CALL mixed_grid_allocate(Gmix,nxGmix,nyGmix) !from agrif_types ! !!!Calculate part of each grid to make Gmix nxG1mix = nxGmix nyG1mix = nyGmix ! !***On rcupre les pts depuis G1 j_min = jmin-1 DO j=1,nyG1mix,(2*rho) i_min = imin-1 DO i=1,nxG1mix,(2*rho) Gmix%nav_lon(i,j) = G1%nav_lon(i_min,j_min) Gmix%nav_lat(i+rho,j) = G1%nav_lat(i_min,j_min) ! Gmix%glam(i,j) = G1%glamt(i_min,j_min) Gmix%glam(i+rho,j) = G1%glamu(i_min,j_min) Gmix%glam(i,j+rho) = G1%glamv(i_min,j_min) Gmix%glam(i+rho,j+rho) = G1%glamf(i_min,j_min) ! Gmix%gphi(i,j) = G1%gphit(i_min,j_min) Gmix%gphi(i+rho,j) = G1%gphiu(i_min,j_min) Gmix%gphi(i,j+rho) = G1%gphiv(i_min,j_min) Gmix%gphi(i+rho,j+rho) = G1%gphif(i_min,j_min) ! Gmix%e1(i,j) = G1%e1t(i_min,j_min) Gmix%e1(i+rho,j) = G1%e1u(i_min,j_min) Gmix%e1(i,j+rho) = G1%e1v(i_min,j_min) Gmix%e1(i+rho,j+rho) = G1%e1f(i_min,j_min) ! Gmix%e2(i,j) = G1%e2t(i_min,j_min) Gmix%e2(i+rho,j) = G1%e2u(i_min,j_min) Gmix%e2(i,j+rho) = G1%e2v(i_min,j_min) Gmix%e2(i+rho,j+rho) = G1%e2f(i_min,j_min) ! i_min = i_min+1 ENDDO j_min = j_min+1 ENDDO ! !print*, 'Gmix%nav_lat= ' !DO j=1,nyG1mix !print*, Gmix%nav_lat(:,j) !END DO ! print*, 'G1%glamt= ' DO j=jmin-1,jmax+1 print*, G1%glamt(imin-1:imax+1,j) ENDDO print*,'' print*, 'G1%glamu= ' DO j=jmin-1,jmax+1 print*, G1%glamu(imin-1:imax+1,j) ENDDO print*,'' print*, 'Gmix%glam= ' DO j=1,nyGmix print*, Gmix%glam(:,j) ENDDO ! END SUBROUTINE ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE mixed_G2(G1,G2,Gmix) ! TYPE(Coordinates), INTENT(IN) :: G1,G2 TYPE(mixed_coordinates), INTENT(INOUT) :: Gmix ! WRITE(*,*) '' WRITE(*,*) '*** ROUTINE mixed_G2 ***' WRITE(*,*) '' ! ! nxG1 = SIZE(G1%glamt,1) - (imin-1) + 1 !(+1) rajout d'une bande fantme nyG1 = (jmax+1) - (jmin-1) + 1 print*, nxG1, nyG1 ! nxG2 = (imax+1) - 2 !(-2) on supprime les bandes de chevauchement nyG2 = nyG1 print*, nxG2, nyG2 ! nx_coarse = nxG1+nxG2 ny_coarse = nyG1 ! !!!Calculate size of mixed grid (nx x ny) nxGmix = (nx_coarse) * 2 !nbre de pts connus (T,U,V,F) suivant x nxGmix = nxGmix + (rho-1)*(nxGmix-1) !nbre de points interpoler ! nyGmix = (ny_coarse) * 2 !nbre de pts connus (T,U,V,F) suivant y nyGmix = nyGmix + (rho-1)*(nyGmix-1) !nbre de points interpoler ! WRITE(*,*) '*** SIZE OF MIXED GRID ***' WRITE(*,*) nxGmix, ' x ', nyGmix WRITE(*,*) '' ! !nx_fine = (nxG1-1 + nxG2-1)*rho + 1 !(-1) on supprime les 2 bandes fantmes !ny_fine = (nyG1-1)*rho + 1 nx_fine = (nx_coarse-2)*rho + 1 ny_fine = (ny_coarse-2)*rho + 1 ! WRITE(*,*) '*** SIZE OF FINE GRID ***' WRITE(*,*) nx_fine, ' x ', ny_fine WRITE(*,*) '' ! !!!Allocate mixed grid CALL mixed_grid_allocate(Gmix,nxGmix,nyGmix) !from agrif_types ! !!!Calculate part of each grid to make Gmix nxG1mix = nxG1 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG1mix = nxG1mix + (rho-1)*(nxG1mix-1) !nbre de points interpoler ! nyG1mix = nyG1 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG1mix = nyG1mix + (rho-1)*(nyG1mix-1) !nbre de points interpoler print*,'nxG1mix= ', nxG1mix, 'nyG1mix= ', nyG1mix ! nxG2mix = nxG2 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG2mix = nxG2mix + (rho-1)*(nxG2mix-1) + 1 !nbre de points interpoler ! nyG2mix = nyG2 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG2mix = nyG2mix + (rho-1)*(nyG2mix-1) !nbre de points interpoler print*,'nxG2mix= ', nxG2mix, 'nyG2mix= ', nyG2mix ! IF((nxG1mix+nxG2mix).NE.nxGmix) THEN WRITE(*,*) '' WRITE(*,*) '*** ERROR ***' WRITE(*,*) 'nxG1mix + nxG2mix /= nxGmix' WRITE(*,*) nxG1mix + nxG2mix, ' /= ', nxGmix WRITE(*,*) '' ENDIF ! !*** On rcupre les pts depuis G1 j_min = jmin-1 DO j=1,nyG1mix,(2*rho) i_min = imin-1 DO i=1,nxG1mix,(2*rho) Gmix%nav_lon(i,j) = G1%nav_lon(i_min,j_min) Gmix%nav_lat(i+rho,j) = G1%nav_lat(i_min,j_min) ! Gmix%glam(i,j) = G1%glamt(i_min,j_min) Gmix%glam(i+rho,j) = G1%glamu(i_min,j_min) Gmix%glam(i,j+rho) = G1%glamv(i_min,j_min) Gmix%glam(i+rho,j+rho) = G1%glamf(i_min,j_min) ! Gmix%gphi(i,j) = G1%gphit(i_min,j_min) Gmix%gphi(i+rho,j) = G1%gphiu(i_min,j_min) Gmix%gphi(i,j+rho) = G1%gphiv(i_min,j_min) Gmix%gphi(i+rho,j+rho) = G1%gphif(i_min,j_min) ! Gmix%e1(i,j) = G1%e1t(i_min,j_min) Gmix%e1(i+rho,j) = G1%e1u(i_min,j_min) Gmix%e1(i,j+rho) = G1%e1v(i_min,j_min) Gmix%e1(i+rho,j+rho) = G1%e1f(i_min,j_min) ! Gmix%e2(i,j) = G1%e2t(i_min,j_min) Gmix%e2(i+rho,j) = G1%e2u(i_min,j_min) Gmix%e2(i,j+rho) = G1%e2v(i_min,j_min) Gmix%e2(i+rho,j+rho) = G1%e2f(i_min,j_min) ! i_min = i_min+1 ENDDO j_min = j_min+1 ENDDO ! !***On rcupre les pts depuis G2 j_min = jmin - 1 DO j=1,nyG2mix,(2*rho) i_min = 3 DO k=i,nxGmix,(2*rho) print*, k Gmix%nav_lon(k,j) = G2%nav_lon(i_min,j_min) Gmix%nav_lat(k+rho,j) = G2%nav_lat(i_min,j_min) ! Gmix%glam(k,j) = G2%glamt(i_min,j_min) Gmix%glam(k+rho,j) = G2%glamu(i_min,j_min) Gmix%glam(k,j+rho) = G2%glamv(i_min,j_min) Gmix%glam(k+rho,j+rho) = G2%glamf(i_min,j_min) ! Gmix%gphi(k,j) = G2%gphit(i_min,j_min) Gmix%gphi(k+rho,j) = G2%gphiu(i_min,j_min) Gmix%gphi(k,j+rho) = G2%gphiv(i_min,j_min) Gmix%gphi(k+rho,j+rho) = G2%gphif(i_min,j_min) ! Gmix%e1(k,j) = G2%e1t(i_min,j_min) Gmix%e1(k+rho,j) = G2%e1u(i_min,j_min) Gmix%e1(k,j+rho) = G2%e1v(i_min,j_min) Gmix%e1(k+rho,j+rho) = G2%e1f(i_min,j_min) ! Gmix%e2(k,j) = G2%e2t(i_min,j_min) Gmix%e2(k+rho,j) = G2%e2u(i_min,j_min) Gmix%e2(k,j+rho) = G2%e2v(i_min,j_min) Gmix%e2(k+rho,j+rho) = G2%e2f(i_min,j_min) ! i_min = i_min+1 ENDDO j_min = j_min+1 ENDDO ! print*, 'G1%glamt= ' DO j=jmin-1,jmax+1 print*, G1%glamt(imin-1:SIZE(G1%glamt,1),j) ENDDO print*,'' print*, 'G1%glamu= ' DO j=jmin-1,jmax+1 print*, G1%glamu(imin-1:SIZE(G1%glamu,1),j) ENDDO print*,'' print*, 'Gmix%glam= ' DO j=1,nyGmix print*, Gmix%glam(:,j) ENDDO END SUBROUTINE ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE mixed_G3(G1,G3,Gmix) ! TYPE(Coordinates), INTENT(IN) :: G1,G2,G3 TYPE(mixed_coordinates), INTENT(INOUT) :: Gmix ! WRITE(*,*) '' WRITE(*,*) '*** ROUTINE mixed_G3 ***' WRITE(*,*) '' ! ! nxG1 = (imax+1) - (imin-1) + 1 !(+1) rajout de 2 bandes fantmes nyG1 = SIZE(G1%glamt,2) - (jmin-1) + 1 print*, nxG1, nyG1 ! nxG3 = nxG1 nyG3 = SIZE(G3%glamt,2) - (jmax-1) + 1 - 2 !(-2) on supprime les bandes de chevauchement print*, nxG3, nyG3 ! nx_coarse = nxG1 ny_coarse = nyG1 + nyG3 ! !!!Calculate size of mixed grid (nx x ny) nxGmix = (nx_coarse) * 2 !nbre de pts connus (T,U,V,F) suivant x nxGmix = nxGmix + (rho-1)*(nxGmix-1) !nbre de points interpoler ! nyGmix = (ny_coarse) * 2 !nbre de pts connus (T,U,V,F) suivant y nyGmix = nyGmix + (rho-1)*(nyGmix-1) !nbre de points interpoler ! WRITE(*,*) '*** SIZE OF MIXED GRID ***' WRITE(*,*) nxGmix, ' x ', nyGmix WRITE(*,*) '' ! !nx_fine = (nxGmix/2) - 1 !ny_fine = (nyGmix/2) - 1 nx_fine = (nx_coarse-2)*rho + 1 ny_fine = (ny_coarse-2)*rho + 1 ! WRITE(*,*) '*** SIZE OF FINE GRID ***' WRITE(*,*) nx_fine, ' x ', ny_fine WRITE(*,*) '' ! !!!Allocate mixed grid CALL mixed_grid_allocate(Gmix,nxGmix,nyGmix) !from agrif_types ! !!!Calculate part of each grid to make Gmix nxG1mix = nxG1 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG1mix = nxG1mix + (rho-1)*(nxG1mix-1) !nbre de points interpoler ! nyG1mix = nyG1 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG1mix = nyG1mix + (rho-1)*(nyG1mix-1) !nbre de points interpoler print*,'nxG1mix= ', nxG1mix, 'nyG1mix= ', nyG1mix ! nxG2mix = nxG2 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG2mix = nxG2mix + (rho-1)*(nxG2mix-1) !nbre de points interpoler ! nyG2mix = nyG2 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG2mix = nyG2mix + (rho-1)*(nyG2mix-1) !nbre de points interpoler print*,'nxG2mix= ', nxG2mix, 'nyG2mix= ', nyG2mix ! nxG3mix = nxG3 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG3mix = nxG3mix + (rho-1)*(nxG3mix-1) !nbre de points interpoler ! nyG3mix = nyG3 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG3mix = nyG3mix + (rho-1)*(nyG3mix-1) !nbre de points interpoler print*,'nxG3mix= ', nxG3mix, 'nyG3mix= ', nyG3mix ! nxG4mix = nxG4 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG4mix = nxG4mix + (rho-1)*(nxG4mix-1) !nbre de points interpoler ! nyG4mix = nyG4 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG4mix = nyG4mix + (rho-1)*(nyG4mix-1) !nbre de points interpoler print*,'nxG4mix= ', nxG4mix, 'nyG4mix= ', nyG4mix ! !***On rcupre les pts depuis G1 j_min = jmin-1 DO j=1,nyG1mix,(2*rho) i_min = imin-1 DO i=1,nxG1mix,(2*rho) Gmix%nav_lon(i,j) = G1%nav_lon(i_min,j_min) Gmix%nav_lat(i+rho,j) = G1%nav_lat(i_min,j_min) ! Gmix%glam(i,j) = G1%glamt(i_min,j_min) Gmix%glam(i+rho,j) = G1%glamu(i_min,j_min) Gmix%glam(i,j+rho) = G1%glamv(i_min,j_min) Gmix%glam(i+rho,j+rho) = G1%glamf(i_min,j_min) ! Gmix%gphi(i,j) = G1%gphit(i_min,j_min) Gmix%gphi(i+rho,j) = G1%gphiu(i_min,j_min) Gmix%gphi(i,j+rho) = G1%gphiv(i_min,j_min) Gmix%gphi(i+rho,j+rho) = G1%gphif(i_min,j_min) ! Gmix%e1(i,j) = G1%e1t(i_min,j_min) Gmix%e1(i+rho,j) = G1%e1u(i_min,j_min) Gmix%e1(i,j+rho) = G1%e1v(i_min,j_min) Gmix%e1(i+rho,j+rho) = G1%e1f(i_min,j_min) ! Gmix%e2(i,j) = G1%e2t(i_min,j_min) Gmix%e2(i+rho,j) = G1%e2u(i_min,j_min) Gmix%e2(i,j+rho) = G1%e2v(i_min,j_min) Gmix%e2(i+rho,j+rho) = G1%e2f(i_min,j_min) ! i_min = i_min+1 ENDDO j_min = j_min+1 ENDDO ! !***On rcupre les pts depuis G2 j_min = jmin - 1 DO j=1,nyG2mix,(2*rho) i_min = 3 DO k=i,nxGmix,(2*rho) Gmix%nav_lon(k,j) = G2%nav_lon(i_min,j_min) Gmix%nav_lat(k+rho,j) = G2%nav_lat(i_min,j_min) ! Gmix%glam(k,j) = G2%glamt(i_min,j_min) Gmix%glam(k+rho,j) = G2%glamu(i_min,j_min) Gmix%glam(k,j+rho) = G2%glamv(i_min,j_min) Gmix%glam(k+rho,j+rho) = G2%glamf(i_min,j_min) ! Gmix%gphi(k,j) = G2%gphit(i_min,j_min) Gmix%gphi(k+rho,j) = G2%gphiu(i_min,j_min) Gmix%gphi(k,j+rho) = G2%gphiv(i_min,j_min) Gmix%gphi(k+rho,j+rho) = G2%gphif(i_min,j_min) ! Gmix%e1(k,j) = G2%e1t(i_min,j_min) Gmix%e1(k+rho,j) = G2%e1u(i_min,j_min) Gmix%e1(k,j+rho) = G2%e1v(i_min,j_min) Gmix%e1(k+rho,j+rho) = G2%e1f(i_min,j_min) ! Gmix%e2(k,j) = G2%e2t(i_min,j_min) Gmix%e2(k+rho,j) = G2%e2u(i_min,j_min) Gmix%e2(k,j+rho) = G2%e2v(i_min,j_min) Gmix%e2(k+rho,j+rho) = G2%e2f(i_min,j_min) ! i_min = i_min+1 ENDDO j_min = j_min+1 ENDDO ! !**On rcupre les pts depuis G3 j_min = SIZE(G3%glamt,2)-2 !on supprime les bandes de chevauchement DO m=j,nyGmix,(2*rho) i_min = imin + imax DO i=1,nxG3mix,(2*rho) Gmix%nav_lon(i,m) = G3%nav_lon(i_min,j_min) Gmix%nav_lat(i+rho,m) = G3%nav_lat(i_min,j_min) ! Gmix%glam(i,m) = G3%glamt(i_min,j_min) Gmix%glam(i+rho,m) = G3%glamu(i_min,j_min) Gmix%glam(i,m+rho) = G3%glamv(i_min,j_min) Gmix%glam(i+rho,m+rho) = G3%glamf(i_min,j_min) ! Gmix%gphi(i,m) = G3%gphit(i_min,j_min) Gmix%gphi(i+rho,m) = G3%gphiu(i_min,j_min) Gmix%gphi(i,m+rho) = G3%gphiv(i_min,j_min) Gmix%gphi(i+rho,m+rho) = G3%gphif(i_min,j_min) ! Gmix%e1(i,m) = G3%e1t(i_min,j_min) Gmix%e1(i+rho,m) = G3%e1u(i_min,j_min) Gmix%e1(i,m+rho) = G3%e1v(i_min,j_min) Gmix%e1(i+rho,m+rho) = G3%e1f(i_min,j_min) ! Gmix%e2(i,m) = G3%e2t(i_min,j_min) Gmix%e2(i+rho,m) = G3%e2u(i_min,j_min) Gmix%e2(i,m+rho) = G3%e2v(i_min,j_min) Gmix%e2(i+rho,m+rho) = G3%e2f(i_min,j_min) ! i_min = i_min - 1 !on se dplace le long de (-i) ENDDO j_min = j_min - 1 !on se dplace le long de (-j) ENDDO ! !**On rcupre les pts depuis G4 j_min = SIZE(G4%glamt,2)-2 !on supprime les bandes de chevauchement DO m=j,nyGmix,(2*rho) i_min = SIZE(G4%glamt,1)-2 DO k=i,nxGmix,(2*rho) Gmix%nav_lon(k,m) = G4%nav_lon(i_min,j_min) Gmix%nav_lat(k+rho,m) = G4%nav_lat(i_min,j_min) ! Gmix%glam(k,m) = G4%glamt(i_min,j_min) Gmix%glam(k+rho,m) = G4%glamu(i_min,j_min) Gmix%glam(k,m+rho) = G4%glamv(i_min,j_min) Gmix%glam(k+rho,m+rho) = G4%glamf(i_min,j_min) ! Gmix%gphi(k,m) = G4%gphit(i_min,j_min) Gmix%gphi(k+rho,m) = G4%gphiu(i_min,j_min) Gmix%gphi(k,m+rho) = G4%gphiv(i_min,j_min) Gmix%gphi(k+rho,m+rho) = G4%gphif(i_min,j_min) ! Gmix%e1(k,m) = G4%e1t(i_min,j_min) Gmix%e1(k+rho,m) = G4%e1u(i_min,j_min) Gmix%e1(k,m+rho) = G4%e1v(i_min,j_min) Gmix%e1(k+rho,m+rho) = G4%e1f(i_min,j_min) ! Gmix%e2(k,m) = G4%e2t(i_min,j_min) Gmix%e2(k+rho,m) = G4%e2u(i_min,j_min) Gmix%e2(k,m+rho) = G4%e2v(i_min,j_min) Gmix%e2(k+rho,m+rho) = G4%e2f(i_min,j_min) ! i_min = i_min - 1 !on se dplace le long de (-i) ENDDO j_min = j_min - 1 !on se dplace le long de (-j) ENDDO ! print*, 'G1%glamt= ' DO j=jmin-1,SIZE(G1%glamt,2) print*, G1%glamt(imin-1:SIZE(G1%glamt,1),j) ENDDO print*,'' print*, 'G1%glamu= ' DO j=jmin-1,SIZE(G1%glamu,2) print*, G1%glamu(imin-1:SIZE(G1%glamu,1),j) ENDDO print*,'' print*, 'Gmix%glam= ' DO j=1,nyGmix print*, Gmix%glam(:,j) ENDDO END SUBROUTINE ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE mixed_G4(G1,G2,G3,G4,Gmix) ! TYPE(Coordinates), INTENT(IN) :: G1,G2,G3,G4 TYPE(mixed_coordinates), INTENT(INOUT) :: Gmix ! WRITE(*,*) '' WRITE(*,*) '*** ROUTINE mixed_G4 ***' WRITE(*,*) '' ! ! nxG1 = SIZE(G1%glamt,1) - (imin-1) + 1 !(+1) rajout d'une bande fantme nyG1 = SIZE(G1%glamt,2) - (jmin-1) + 1 print*, nxG1, nyG1 ! nxG2 = imax + 1 - 2 !(-2) on supprime les bandes de chevauchement nyG2 = SIZE(G2%glamt,2) - (jmin-1) + 1 print*, nxG2, nyG2 ! nxG3 = SIZE(G3%glamt,1) - (imin-1) + 1 nyG3 = SIZE(G3%glamt,2) - (jmax-1) + 1 - 2 !(-2) on supprime les bandes de chevauchement print*, nxG3, nyG3 ! nxG4 = imax + 1 - 2 !(-2) on supprime les bandes de chevauchement nyG4 = SIZE(G4%glamt,2) - (jmax-1) + 1 - 2 !(-2) on supprime les bandes de chevauchement print*, nxG4, nyG4 ! WRITE(*,*) '' WRITE(*,*) '*** CHECKING SIZE OF THE MIXED GRID***' WRITE(*,*) ' nxG1 + nxG2 = nxG3 + nxG4 = nx_coarse' WRITE(*,*) nxG1 + nxG2, ' = ', nxG3 + nxG4 WRITE(*,*) '' WRITE(*,*) ' nyG1 + nyG3 = nyG2 + nyG4 = ny_coarse' WRITE(*,*) nyG1 + nyG3, ' = ', nyG2 + nyG4 WRITE(*,*) '' ! nx_coarse = nxG1 + nxG2 ny_coarse = nyG1 + nyG3 ! !!!Calculate size of mixed grid (nx x ny) nxGmix = (nxG1 + nxG2) * 2 !nbre de pts connus (T,U,V,F) suivant x nxGmix = nxGmix + (rho-1)*(nxGmix-1) !nbre de points interpoler ! nyGmix = (nyG1 + nyG3) * 2 !nbre de pts connus (T,U,V,F) suivant y nyGmix = nyGmix + (rho-1)*(nyGmix-1) !nbre de points interpoler ! WRITE(*,*) '*** SIZE OF MIXED GRID ***' WRITE(*,*) nxGmix, ' x ', nyGmix WRITE(*,*) '' ! !nx_fine = (nxGmix/2) - 1 !ny_fine = (nyGmix/2) - 1 nx_fine = (nx_coarse-2)*rho + 1 ny_fine = (ny_coarse-2)*rho + 1 ! WRITE(*,*) '*** SIZE OF FINE GRID ***' WRITE(*,*) nx_fine, ' x ', ny_fine WRITE(*,*) '' ! !!!Allocate mixed grid CALL mixed_grid_allocate(Gmix,nxGmix,nyGmix) !from agrif_types ! !!!Calculate part of each grid to make Gmix nxG1mix = nxG1 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG1mix = nxG1mix + (rho-1)*(nxG1mix-1) !nbre de points interpoler ! nyG1mix = nyG1 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG1mix = nyG1mix + (rho-1)*(nyG1mix-1) !nbre de points interpoler print*,'nxG1mix= ', nxG1mix, 'nyG1mix= ', nyG1mix ! nxG2mix = nxG2 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG2mix = nxG2mix + (rho-1)*(nxG2mix-1) !nbre de points interpoler ! nyG2mix = nyG2 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG2mix = nyG2mix + (rho-1)*(nyG2mix-1) !nbre de points interpoler print*,'nxG2mix= ', nxG2mix, 'nyG2mix= ', nyG2mix ! nxG3mix = nxG3 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG3mix = nxG3mix + (rho-1)*(nxG3mix-1) !nbre de points interpoler ! nyG3mix = nyG3 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG3mix = nyG3mix + (rho-1)*(nyG3mix-1) !nbre de points interpoler print*,'nxG3mix= ', nxG3mix, 'nyG3mix= ', nyG3mix ! nxG4mix = nxG4 * 2 !nbre de pts connus (T,U,V,F) suivant x nxG4mix = nxG4mix + (rho-1)*(nxG4mix-1) !nbre de points interpoler ! nyG4mix = nyG4 * 2 !nbre de pts connus (T,U,V,F) suivant x nyG4mix = nyG4mix + (rho-1)*(nyG4mix-1) !nbre de points interpoler print*,'nxG4mix= ', nxG4mix, 'nyG4mix= ', nyG4mix ! !***On rcupre les pts depuis G1 j_min = jmin-1 DO j=1,nyG1mix,(2*rho) i_min = imin-1 DO i=1,nxG1mix,(2*rho) Gmix%nav_lon(i,j) = G1%nav_lon(i_min,j_min) Gmix%nav_lat(i+rho,j) = G1%nav_lat(i_min,j_min) ! Gmix%glam(i,j) = G1%glamt(i_min,j_min) Gmix%glam(i+rho,j) = G1%glamu(i_min,j_min) Gmix%glam(i,j+rho) = G1%glamv(i_min,j_min) Gmix%glam(i+rho,j+rho) = G1%glamf(i_min,j_min) ! Gmix%gphi(i,j) = G1%gphit(i_min,j_min) Gmix%gphi(i+rho,j) = G1%gphiu(i_min,j_min) Gmix%gphi(i,j+rho) = G1%gphiv(i_min,j_min) Gmix%gphi(i+rho,j+rho) = G1%gphif(i_min,j_min) ! Gmix%e1(i,j) = G1%e1t(i_min,j_min) Gmix%e1(i+rho,j) = G1%e1u(i_min,j_min) Gmix%e1(i,j+rho) = G1%e1v(i_min,j_min) Gmix%e1(i+rho,j+rho) = G1%e1f(i_min,j_min) ! Gmix%e2(i,j) = G1%e2t(i_min,j_min) Gmix%e2(i+rho,j) = G1%e2u(i_min,j_min) Gmix%e2(i,j+rho) = G1%e2v(i_min,j_min) Gmix%e2(i+rho,j+rho) = G1%e2f(i_min,j_min) ! i_min = i_min+1 ENDDO j_min = j_min+1 ENDDO ! !***On rcupre les pts depuis G2 j_min = jmin - 1 DO j=1,nyG2mix,(2*rho) i_min = 3 DO k=i,nxGmix,(2*rho) Gmix%nav_lon(k,j) = G2%nav_lon(i_min,j_min) Gmix%nav_lat(k+rho,j) = G2%nav_lat(i_min,j_min) ! Gmix%glam(k,j) = G2%glamt(i_min,j_min) Gmix%glam(k+rho,j) = G2%glamu(i_min,j_min) Gmix%glam(k,j+rho) = G2%glamv(i_min,j_min) Gmix%glam(k+rho,j+rho) = G2%glamf(i_min,j_min) ! Gmix%gphi(k,j) = G2%gphit(i_min,j_min) Gmix%gphi(k+rho,j) = G2%gphiu(i_min,j_min) Gmix%gphi(k,j+rho) = G2%gphiv(i_min,j_min) Gmix%gphi(k+rho,j+rho) = G2%gphif(i_min,j_min) ! Gmix%e1(k,j) = G2%e1t(i_min,j_min) Gmix%e1(k+rho,j) = G2%e1u(i_min,j_min) Gmix%e1(k,j+rho) = G2%e1v(i_min,j_min) Gmix%e1(k+rho,j+rho) = G2%e1f(i_min,j_min) ! Gmix%e2(k,j) = G2%e2t(i_min,j_min) Gmix%e2(k+rho,j) = G2%e2u(i_min,j_min) Gmix%e2(k,j+rho) = G2%e2v(i_min,j_min) Gmix%e2(k+rho,j+rho) = G2%e2f(i_min,j_min) ! i_min = i_min+1 ENDDO j_min = j_min+1 ENDDO ! !**On rcupre les pts depuis G3 j_min = SIZE(G3%glamt,2)-2 !on supprime les bandes de chevauchement DO m=j,nyGmix,(2*rho) i_min = imin + imax DO i=1,nxG3mix,(2*rho) Gmix%nav_lon(i,m) = G3%nav_lon(i_min,j_min) Gmix%nav_lat(i+rho,m) = G3%nav_lat(i_min,j_min) ! Gmix%glam(i,m) = G3%glamt(i_min,j_min) Gmix%glam(i+rho,m) = G3%glamu(i_min,j_min) Gmix%glam(i,m+rho) = G3%glamv(i_min,j_min) Gmix%glam(i+rho,m+rho) = G3%glamf(i_min,j_min) ! Gmix%gphi(i,m) = G3%gphit(i_min,j_min) Gmix%gphi(i+rho,m) = G3%gphiu(i_min,j_min) Gmix%gphi(i,m+rho) = G3%gphiv(i_min,j_min) Gmix%gphi(i+rho,m+rho) = G3%gphif(i_min,j_min) ! Gmix%e1(i,m) = G3%e1t(i_min,j_min) Gmix%e1(i+rho,m) = G3%e1u(i_min,j_min) Gmix%e1(i,m+rho) = G3%e1v(i_min,j_min) Gmix%e1(i+rho,m+rho) = G3%e1f(i_min,j_min) ! Gmix%e2(i,m) = G3%e2t(i_min,j_min) Gmix%e2(i+rho,m) = G3%e2u(i_min,j_min) Gmix%e2(i,m+rho) = G3%e2v(i_min,j_min) Gmix%e2(i+rho,m+rho) = G3%e2f(i_min,j_min) ! i_min = i_min - 1 !on se dplace le long de (-i) ENDDO j_min = j_min - 1 !on se dplace le long de (-j) ENDDO ! !**On rcupre les pts depuis G4 j_min = SIZE(G4%glamt,2)-2 !on supprime les bandes de chevauchement DO m=j,nyGmix,(2*rho) i_min = SIZE(G4%glamt,1)-2 DO k=i,nxGmix,(2*rho) Gmix%nav_lon(k,m) = G4%nav_lon(i_min,j_min) Gmix%nav_lat(k+rho,m) = G4%nav_lat(i_min,j_min) ! Gmix%glam(k,m) = G4%glamt(i_min,j_min) Gmix%glam(k+rho,m) = G4%glamu(i_min,j_min) Gmix%glam(k,m+rho) = G4%glamv(i_min,j_min) Gmix%glam(k+rho,m+rho) = G4%glamf(i_min,j_min) ! Gmix%gphi(k,m) = G4%gphit(i_min,j_min) Gmix%gphi(k+rho,m) = G4%gphiu(i_min,j_min) Gmix%gphi(k,m+rho) = G4%gphiv(i_min,j_min) Gmix%gphi(k+rho,m+rho) = G4%gphif(i_min,j_min) ! Gmix%e1(k,m) = G4%e1t(i_min,j_min) Gmix%e1(k+rho,m) = G4%e1u(i_min,j_min) Gmix%e1(k,m+rho) = G4%e1v(i_min,j_min) Gmix%e1(k+rho,m+rho) = G4%e1f(i_min,j_min) ! Gmix%e2(k,m) = G4%e2t(i_min,j_min) Gmix%e2(k+rho,m) = G4%e2u(i_min,j_min) Gmix%e2(k,m+rho) = G4%e2v(i_min,j_min) Gmix%e2(k+rho,m+rho) = G4%e2f(i_min,j_min) ! i_min = i_min - 1 !on se dplace le long de (-i) ENDDO j_min = j_min - 1 !on se dplace le long de (-j) ENDDO ! print*, 'G1%glamt= ' DO j=jmin-1,SIZE(G1%glamt,2) print*, G1%glamt(imin-1:SIZE(G1%glamt,1),j) ENDDO print*,'' print*, 'G1%glamu= ' DO j=jmin-1,SIZE(G1%glamu,2) print*, G1%glamu(imin-1:SIZE(G1%glamu,1),j) ENDDO print*,'' print*, 'Gmix%glam= ' DO j=1,nyGmix print*, Gmix%glam(:,j) ENDDO END SUBROUTINE ! ! !******************************************************** ! SUBROUTINE interp * ! * ! calculate polynomial interpolation at 4th order * ! * ! CALLED from create_coordinates.f90 * !******************************************************** ! SUBROUTINE interp_grid(Gmix) ! TYPE(mixed_coordinates) :: Gmix !TYPE(Coordinates) :: Grid1 REAL*8, DIMENSION(rho-1,4) :: S !Array to store the coefficients of Lagrange INTEGER :: i, j, k, status ! WRITE(*,*) '' WRITE(*,*) '*** ROUTINE interp ***' WRITE(*,*) '' ! !!!Calculate coefficients status = pol_coef(S) ! !IF(rho>1) THEN !!!Interpolation along longitude do j=1,SIZE(Gmix%glam,2),rho do i=1,SIZE(Gmix%glam,1),rho do k=1,rho-1 Gmix%glam(i+rho+k,j) = S(k,1) * Gmix%glam(i,j) & + S(k,2) * Gmix%glam(i+1*rho,j) & + S(k,3) * Gmix%glam(i+2*rho,j) & + S(k,4) * Gmix%glam(i+3*rho,j) ! Gmix%gphi(i+rho+k,j) = S(k,1) * Gmix%gphi(i,j) & + S(k,2) * Gmix%gphi(i+1*rho,j) & + S(k,3) * Gmix%gphi(i+2*rho,j) & + S(k,4) * Gmix%gphi(i+3*rho,j) ! Gmix%e1(i+rho+k,j) = S(k,1) * Gmix%e1(i,j) & + S(k,2) * Gmix%e1(i+1*rho,j) & + S(k,3) * Gmix%e1(i+2*rho,j) & + S(k,4) * Gmix%e1(i+3*rho,j) ! Gmix%e2(i+rho+k,j) = S(k,1) * Gmix%e2(i,j) & + S(k,2) * Gmix%e2(i+1*rho,j) & + S(k,3) * Gmix%e2(i+2*rho,j) & + S(k,4) * Gmix%e2(i+3*rho,j) ! Gmix%nav_lon(i+rho+k,j) = Gmix%glam(i+rho+k,j) Gmix%nav_lat(i+rho+k,j) = Gmix%gphi(i+rho+k,j) ! end do ! IF(i+3*rho == SIZE(Gmix%glam,1)) EXIT end do end do ! !!!Interpolation along latitude do i=1,SIZE(Gmix%glam,1) do j=1,SIZE(Gmix%glam,2),rho do k=1,rho-1 Gmix%glam(i,j+rho+k) = S(k,1) * Gmix%glam(i,j) & + S(k,2) * Gmix%glam(i,j+rho) & + S(k,3) * Gmix%glam(i,j+2*rho) & + S(k,4) * Gmix%glam(i,j+3*rho) ! Gmix%gphi(i,j+rho+k) = S(k,1) * Gmix%gphi(i,j) & + S(k,2) * Gmix%gphi(i,j+rho) & + S(k,3) * Gmix%gphi(i,j+2*rho) & + S(k,4) * Gmix%gphi(i,j+3*rho) ! Gmix%e1(i,j+rho+k) = S(k,1) * Gmix%e1(i,j) & + S(k,2) * Gmix%e1(i,j+rho) & + S(k,3) * Gmix%e1(i,j+2*rho) & + S(k,4) * Gmix%e1(i,j+3*rho) ! Gmix%e2(i,j+rho+k) = S(k,1) * Gmix%e2(i,j) & + S(k,2) * Gmix%e2(i,j+rho) & + S(k,3) * Gmix%e2(i,j+2*rho) & + S(k,4) * Gmix%e2(i,j+3*rho) ! Gmix%nav_lon(i,j+rho+k) = Gmix%glam(i,j+rho+k) Gmix%nav_lat(i,j+rho+k) = Gmix%gphi(i,j+rho+k) ! end do ! IF(j+3*rho == SIZE(Gmix%glam,2)) EXIT end do end do ! END SUBROUTINE ! ! !******************************************************** ! FUNCTION pol_coef * ! * ! calculate the coefficients of Lagrange * ! for the polynomial interpolation * ! * ! CALLED from SUBROUTINE interp * !******************************************************** ! REAL FUNCTION pol_coef(vect) ! REAL*8, DIMENSION(rho-1,4) :: vect REAL*8, DIMENSION(3) :: v INTEGER :: i, m, k REAL*8 :: x0 !position relative du point calculer INTEGER :: x_k, x_i !position relative des points utiliss pour l'interpolation REAL*8 :: eps ! !on parle de position relative puisque nous utilisons les positions !indiciaires, lesquelles sont rptes dans toute la grille. !Il n'est donc ncessaire de calculer qu'une fois les 4 coefficients !qui seront utiliss dans toute la grille en fonction de rho ! ! WRITE(*,*) '' WRITE(*,*) '*** FUNCTION pol_coef ***' WRITE(*,*) '' ! m=1 eps = 1.-1e-8 ! do k=1,rho-1 x0=rho+1+k i=1 do x_i=1,4+3*(rho-1),rho m=1 do x_k=1,4+3*(rho-1),rho IF(x_k == x_i) THEN CYCLE ELSE v(m) = (x0-x_k) / (x_i-x_k) m=m+1 END IF end do vect(k,i) = product(v) i=i+1 end do end do ! IF((SUM(vect)eps) .AND. rho>1) THEN print*,'' print*, '*** CHECK LAGRANGE COEFFICIENTS: ***' print*,'' do i=1,rho-1 print*,'point #',i print*, 'S(i,:)= ', vect(i,:) print*, 'SUM(S(i,:)) =', SUM(vect(i,:)) print*,'' end do print*, '' pol_coef = 1 ELSE !print*,'Error with Lagrange coefficients' pol_coef = 0 END IF ! END FUNCTION pol_coef ! ! !******************************************************** ! SUBROUTINE alloc_child_grid * ! * ! create the child grid from mixed grid * ! * ! CALLED from create_coordinates.f90 * !******************************************************** ! SUBROUTINE alloc_child_grid(Gmix, Grid1) ! TYPE(mixed_coordinates), INTENT(IN) :: Gmix TYPE(coordinates), INTENT(INOUT) :: Grid1 INTEGER :: i, j, p, q ! WRITE(*,*) '' WRITE(*,*) '*** ROUTINE alloc_child_grid ***' WRITE(*,*) '' ! ! IF(rho>1) THEN q=1 DO j=rho+2,(SIZE(Gmix%glam,2)-(rho)),2 p=1 DO i=rho+2,(SIZE(Gmix%glam,1)-(rho)),2 ! Grid1%nav_lon(p,q) = Gmix%nav_lon(i,j) Grid1%nav_lat(p,q) = Gmix%nav_lat(i,j) ! Grid1%glamt(p,q) = Gmix%glam(i,j) Grid1%glamu(p,q) = Gmix%glam(i+1,j) Grid1%glamv(p,q) = Gmix%glam(i,j+1) Grid1%glamf(p,q) = Gmix%glam(i+1,j+1) ! Grid1%gphit(p,q) = Gmix%gphi(i,j) Grid1%gphiu(p,q) = Gmix%gphi(i+1,j) Grid1%gphiv(p,q) = Gmix%gphi(i,j+1) Grid1%gphif(p,q) = Gmix%gphi(i+1,j+1) ! Grid1%e1t(p,q) = Gmix%e1(i,j) Grid1%e1u(p,q) = Gmix%e1(i+1,j) Grid1%e1v(p,q) = Gmix%e1(i,j+1) Grid1%e1f(p,q) = Gmix%e1(i+1,j+1) ! Grid1%e2t(p,q) = Gmix%e2(i,j) Grid1%e2u(p,q) = Gmix%e2(i+1,j) Grid1%e2v(p,q) = Gmix%e2(i,j+1) Grid1%e2f(p,q) = Gmix%e2(i+1,j+1) ! p=p+1 END DO q=q+1 END DO ELSE q=1 DO j=1,SIZE(Gmix%glam,2),2 p=1 DO i=1,SIZE(Gmix%glam,1),2 ! Grid1%nav_lon(p,q) = Gmix%nav_lon(i,j) Grid1%nav_lat(p,q) = Gmix%nav_lat(i+1,j) ! Grid1%glamt(p,q) = Gmix%glam(i,j) Grid1%glamu(p,q) = Gmix%glam(i+1,j) Grid1%glamv(p,q) = Gmix%glam(i,j+1) Grid1%glamf(p,q) = Gmix%glam(i+1,j+1) ! Grid1%gphit(p,q) = Gmix%gphi(i,j) Grid1%gphiu(p,q) = Gmix%gphi(i+1,j) Grid1%gphiv(p,q) = Gmix%gphi(i,j+1) Grid1%gphif(p,q) = Gmix%gphi(i+1,j+1) ! Grid1%e1t(p,q) = Gmix%e1(i,j) Grid1%e1u(p,q) = Gmix%e1(i+1,j) Grid1%e1v(p,q) = Gmix%e1(i,j+1) Grid1%e1f(p,q) = Gmix%e1(i+1,j+1) ! Grid1%e2t(p,q) = Gmix%e2(i,j) Grid1%e2u(p,q) = Gmix%e2(i+1,j) Grid1%e2v(p,q) = Gmix%e2(i,j+1) Grid1%e2f(p,q) = Gmix%e2(i+1,j+1) ! p=p+1 END DO q=q+1 END DO ENDIF WRITE(*,*) '' WRITE(*,*) 'Size of the Child grid: ', nx_fine, ' x ', ny_fine WRITE(*,*) '' print*, 'Grid1%nav_lat= ' DO j=1,ny_fine print*, Grid1%nav_lat(:,j) END DO ! END SUBROUTINE alloc_child_grid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE tools_brice