SUBROUTINE ana_grid (ng, tile, model) ! !! git $Id$ !! svn $Id: ana_grid.h 1151 2023-02-09 03:08:53Z arango $ !!====================================================================== !! Copyright (c) 2002-2023 The ROMS/TOMS Group ! !! Licensed under a MIT/X style license ! !! See License_ROMS.md ! !======================================================================= ! ! ! This routine sets model grid using an analytical expressions. ! ! ! ! On Output: stored in common blocks: ! ! ! ! "grid" (file grid.h) ! ! "scalars" (file scalar.h) ! ! ! ! el Length (m) of domain box in the ETA-direction. ! ! f Coriolis parameter (1/seconds) at RHO-points. ! ! h Bathymetry (meters; positive) at RHO-points. ! ! hmin Minimum depth of bathymetry (m). ! ! hmax Maximum depth of bathymetry (m). ! ! pm Coordinate transformation metric "m" (1/meters) ! ! associated with the differential distances in XI ! ! at RHO-points. ! ! pn Coordinate transformation metric "n" (1/meters) ! ! associated with the differential distances in ETA. ! ! at RHO-points. ! ! xl Length (m) of domain box in the XI-direction. ! ! xp XI-coordinates (m) at PSI-points. ! ! xr XI-coordinates (m) at RHO-points. ! ! yp ETA-coordinates (m) at PSI-points. ! ! yr ETA-coordinates (m) at RHO-points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_ncparam ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__ ! #include "tile.h" ! CALL ana_grid_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & GRID(ng) % angler, & #if defined CURVGRID && defined UV_ADV & GRID(ng) % dmde, & & GRID(ng) % dndx, & #endif #ifdef ICESHELF & GRID(ng) % zice, & #endif #ifdef SPHERICAL & GRID(ng) % lonp, & & GRID(ng) % lonr, & & GRID(ng) % lonu, & & GRID(ng) % lonv, & & GRID(ng) % latp, & & GRID(ng) % latr, & & GRID(ng) % latu, & & GRID(ng) % latv, & #else & GRID(ng) % xp, & & GRID(ng) % xr, & & GRID(ng) % xu, & & GRID(ng) % xv, & & GRID(ng) % yp, & & GRID(ng) % yr, & & GRID(ng) % yu, & & GRID(ng) % yv, & #endif & GRID(ng) % pn, & & GRID(ng) % pm, & & GRID(ng) % f, & & GRID(ng) % h) ! ! Set analytical header file name used. ! #ifdef DISTRIBUTE IF (Lanafile) THEN #else IF (Lanafile.and.(tile.eq.0)) THEN #endif ANANAME( 7)=MyFile END IF ! RETURN END SUBROUTINE ana_grid ! !*********************************************************************** SUBROUTINE ana_grid_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & angler, & #if defined CURVGRID && defined UV_ADV & dmde, dndx, & #endif #ifdef ICESHELF & zice, & #endif #ifdef SPHERICAL & lonp, lonr, lonu, lonv, & & latp, latr, latu, latv, & #else & xp, xr, xu, xv, & & yp, yr, yu, yv, & #endif & pn, pm, f, h) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_scalars ! #ifdef DISTRIBUTE USE distribute_mod, ONLY : mp_reduce #endif USE exchange_2d_mod, ONLY : exchange_r2d_tile #ifdef DISTRIBUTE USE mp_exchange_mod, ONLY : mp_exchange2d #endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS ! #ifdef ASSUMED_SHAPE real(r8), intent(out) :: angler(LBi:,LBj:) # if defined CURVGRID && defined UV_ADV real(r8), intent(out) :: dmde(LBi:,LBj:) real(r8), intent(out) :: dndx(LBi:,LBj:) # endif # ifdef ICESHELF real(r8), intent(out) :: zice(LBi:,LBj:) # endif # ifdef SPHERICAL real(r8), intent(out) :: lonp(LBi:,LBj:) real(r8), intent(out) :: lonr(LBi:,LBj:) real(r8), intent(out) :: lonu(LBi:,LBj:) real(r8), intent(out) :: lonv(LBi:,LBj:) real(r8), intent(out) :: latp(LBi:,LBj:) real(r8), intent(out) :: latr(LBi:,LBj:) real(r8), intent(out) :: latu(LBi:,LBj:) real(r8), intent(out) :: latv(LBi:,LBj:) # else real(r8), intent(out) :: xp(LBi:,LBj:) real(r8), intent(out) :: xr(LBi:,LBj:) real(r8), intent(out) :: xu(LBi:,LBj:) real(r8), intent(out) :: xv(LBi:,LBj:) real(r8), intent(out) :: yp(LBi:,LBj:) real(r8), intent(out) :: yr(LBi:,LBj:) real(r8), intent(out) :: yu(LBi:,LBj:) real(r8), intent(out) :: yv(LBi:,LBj:) # endif real(r8), intent(out) :: pn(LBi:,LBj:) real(r8), intent(out) :: pm(LBi:,LBj:) real(r8), intent(out) :: f(LBi:,LBj:) real(r8), intent(out) :: h(LBi:,LBj:) #else real(r8), intent(out) :: angler(LBi:UBi,LBj:UBj) # if defined CURVGRID && defined UV_ADV real(r8), intent(out) :: dmde(LBi:UBi,LBj:UBj) real(r8), intent(out) :: dndx(LBi:UBi,LBj:UBj) # endif # ifdef ICESHELF real(r8), intent(out) :: zice(LBi:UBi,LBj:UBj) # endif # ifdef SPHERICAL real(r8), intent(out) :: lonp(LBi:UBi,LBj:UBj) real(r8), intent(out) :: lonr(LBi:UBi,LBj:UBj) real(r8), intent(out) :: lonu(LBi:UBi,LBj:UBj) real(r8), intent(out) :: lonv(LBi:UBi,LBj:UBj) real(r8), intent(out) :: latp(LBi:UBi,LBj:UBj) real(r8), intent(out) :: latr(LBi:UBi,LBj:UBj) real(r8), intent(out) :: latu(LBi:UBi,LBj:UBj) real(r8), intent(out) :: latv(LBi:UBi,LBj:UBj) # else real(r8), intent(out) :: xp(LBi:UBi,LBj:UBj) real(r8), intent(out) :: xr(LBi:UBi,LBj:UBj) real(r8), intent(out) :: xu(LBi:UBi,LBj:UBj) real(r8), intent(out) :: xv(LBi:UBi,LBj:UBj) real(r8), intent(out) :: yp(LBi:UBi,LBj:UBj) real(r8), intent(out) :: yr(LBi:UBi,LBj:UBj) real(r8), intent(out) :: yu(LBi:UBi,LBj:UBj) real(r8), intent(out) :: yv(LBi:UBi,LBj:UBj) # endif real(r8), intent(out) :: pn(LBi:UBi,LBj:UBj) real(r8), intent(out) :: pm(LBi:UBi,LBj:UBj) real(r8), intent(out) :: f(LBi:UBi,LBj:UBj) real(r8), intent(out) :: h(LBi:UBi,LBj:UBj) #endif ! ! Local variable declarations. ! integer :: Imin, Imax, Jmin, Jmax integer :: NSUB, i, j, k ! real(r8) :: Esize, Xsize, beta, depth real(r8) :: dx, dy, f0, my_min, my_max, val1 #ifdef DISTRIBUTE real(r8), dimension(2) :: buffer character (len=3), dimension(2) :: op_handle #endif real(r8) :: wrkX(IminS:ImaxS,JminS:JmaxS) real(r8) :: wrkY(IminS:ImaxS,JminS:JmaxS) #include "set_bounds.h" ! !----------------------------------------------------------------------- ! Set grid parameters: ! ! Xsize Length (m) of domain box in the XI-direction. ! Esize Length (m) of domain box in the ETA-direction. ! depth Maximum depth of bathymetry (m). ! f0 Coriolis parameter, f-plane constant (1/s). ! beta Coriolis parameter, beta-plane constant (1/s/m). !----------------------------------------------------------------------- ! #if defined MY_APPLICATION Xsize=??? Esize=??? depth=??? f0=??? beta=??? #else ana_grid: no values provided for Xsize, Esize, depth, f0, and beta. #endif ! ! Load grid parameters to global storage. ! IF (DOMAIN(ng)%NorthEast_Test(tile)) THEN xl(ng)=Xsize el(ng)=Esize END IF ! !----------------------------------------------------------------------- ! Compute the (XI,ETA) coordinates at PSI- and RHO-points. ! Set grid spacing (m). !----------------------------------------------------------------------- ! ! Determine I- and J-ranges for computing grid data. These ranges ! are special in periodic boundary conditons since periodicity cannot ! be imposed in the grid coordinates. ! IF (DOMAIN(ng)%Western_Edge(tile)) THEN Imin=Istr-1 ELSE Imin=Istr END IF IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN Imax=Iend+1 ELSE Imax=Iend END IF IF (DOMAIN(ng)%Southern_Edge(tile)) THEN Jmin=Jstr-1 ELSE Jmin=Jstr END IF IF (DOMAIN(ng)%Northern_Edge(tile)) THEN Jmax=Jend+1 ELSE Jmax=Jend END IF dx=Xsize/REAL(Lm(ng),r8) dy=Esize/REAL(Mm(ng),r8) DO j=Jmin,Jmax DO i=Imin,Imax xp(i,j)=dx*REAL(i-1,r8) xr(i,j)=dx*(REAL(i-1,r8)+0.5_r8) xu(i,j)=xp(i,j) xv(i,j)=xr(i,j) yp(i,j)=dy*REAL(j-1,r8) yr(i,j)=dy*(REAL(j-1,r8)+0.5_r8) yu(i,j)=yr(i,j) yv(i,j)=yp(i,j) END DO END DO #ifdef DISTRIBUTE ! ! Exchange boundary data. ! CALL mp_exchange2d (ng, tile, model, 4, & & LBi, UBi, LBj, UBj, & & NghostPoints, .FALSE., .FALSE., & & xp, xr, xu, xv) CALL mp_exchange2d (ng, tile, model, 4, & & LBi, UBi, LBj, UBj, & & NghostPoints, .FALSE., .FALSE., & & yp, yr, yu, yv) #endif ! !----------------------------------------------------------------------- ! Compute coordinate transformation metrics at RHO-points "pm" and ! "pn" (1/m) associated with the differential distances in XI and ! ETA, respectively. !----------------------------------------------------------------------- ! #define J_RANGE MIN(JstrT,Jstr-1),MAX(Jend+1,JendT) #define I_RANGE MIN(IstrT,Istr-1),MAX(Iend+1,IendT) DO j=J_RANGE DO i=I_RANGE wrkX(i,j)=1.0_r8/dx wrkY(i,j)=1.0_r8/dy END DO END DO #undef J_RANGE #undef I_RANGE DO j=JstrT,JendT DO i=IstrT,IendT pm(i,j)=wrkX(i,j) pn(i,j)=wrkY(i,j) END DO END DO ! ! Exchange boundary data. ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & pm) CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & pn) END IF #ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, model, 2, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & pm, pn) #endif #if (defined CURVGRID && defined UV_ADV) ! !----------------------------------------------------------------------- ! Compute d(1/n)/d(xi) and d(1/m)/d(eta) at RHO-points. !----------------------------------------------------------------------- ! DO j=Jstr,Jend DO i=Istr,Iend dndx(i,j)=0.5_r8*((1.0_r8/wrkY(i+1,j ))- & & (1.0_r8/wrkY(i-1,j ))) dmde(i,j)=0.5_r8*((1.0_r8/wrkX(i ,j+1))- & & (1.0_r8/wrkX(i ,j-1))) END DO END DO ! ! Exchange boundary data. ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & dndx) CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & dmde) END IF # ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, model, 2, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & dndx, dmde) # endif #endif ! !----------------------------------------------------------------------- ! Angle (radians) between XI-axis and true EAST at RHO-points. !----------------------------------------------------------------------- ! #if defined MY_APPLICATION DO j=JstrT,JendT DO i=IstrT,IendT angler(i,j)=??? END DO END DO #else DO j=JstrT,JendT DO i=IstrT,IendT angler(i,j)=0.0_r8 END DO END DO #endif ! ! Exchange boundary data. ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & angler) END IF #ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, model, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & angler) #endif ! !----------------------------------------------------------------------- ! Compute Coriolis parameter (1/s) at RHO-points. !----------------------------------------------------------------------- ! #if defined MY_APPLICATION DO j=JstrT,JendT DO i=IstrT,IendT f(i,j)=??? END DO END DO #else val1=0.5_r8*Esize DO j=JstrT,JendT DO i=IstrT,IendT f(i,j)=f0+beta*(yr(i,j)-val1) END DO END DO #endif ! ! Exchange boundary data. ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & f) END IF #ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, model, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & f) #endif ! !----------------------------------------------------------------------- ! Set bathymetry (meters; positive) at RHO-points. !----------------------------------------------------------------------- ! #if defined MY_APPLICATION DO j=JstrT,JendT DO i=IstrT,IendT h(i,j)=??? END DO END DO #else DO j=JstrT,JendT DO i=IstrT,IendT h(i,j)=depth END DO END DO #endif ! ! Exchange boundary data. ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & h) END IF #ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, model, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & h) #endif #ifdef ICESHELF ! !----------------------------------------------------------------------- ! Set depth of ice shelf (meters; negative) at RHO-points. !----------------------------------------------------------------------- ! # ifdef WEDDELL DO j=JstrT,JendT DO i=IstrT,IendT zice(i,j)=??? END DO END DO # else ana_grid.h: No value provided for zice. # endif ! ! Exchange boundary data. ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & zice) END IF # ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, model, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & zice) # endif #endif ! ! Determine minimum depth: first, determine minimum values of depth ! within each subdomain, then determine global minimum by comparing ! these subdomain minima. ! my_min=h(IstrT,JstrT) my_max=h(IstrT,JstrT) DO j=JstrT,JendT DO i=IstrT,IendT my_min=MIN(my_min,h(i,j)) my_max=MAX(my_max,h(i,j)) END DO END DO #ifdef DISTRIBUTE NSUB=1 ! distributed-memory #else IF (DOMAIN(ng)%SouthWest_Corner(tile).and. & & DOMAIN(ng)%NorthEast_Corner(tile)) THEN NSUB=1 ! non-tiled application ELSE NSUB=NtileX(ng)*NtileE(ng) ! tiled application END IF #endif !$OMP CRITICAL (H_RANGE) IF (tile_count.eq.0) THEN hmin(ng)=my_min hmax(ng)=my_max ELSE hmin(ng)=MIN(hmin(ng),my_min) hmax(ng)=MAX(hmax(ng),my_max) END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 #ifdef DISTRIBUTE buffer(1)=hmin(ng) buffer(2)=hmax(ng) op_handle(1)='MIN' op_handle(2)='MAX' CALL mp_reduce (ng, model, 2, buffer, op_handle) hmin(ng)=buffer(1) hmax(ng)=buffer(2) #endif END IF !$OMP END CRITICAL (H_RANGE) ! RETURN END SUBROUTINE ana_grid_tile