MODULE metrics_mod ! !svn $Id: metrics.F 352 2009-05-29 20:57:39Z arango $ !======================================================================= ! Copyright (c) 2002-2009 The ROMS/TOMS Group ! ! Licensed under a MIT/X style license ! ! See License_ROMS.txt Hernan G. Arango ! !========================================== Alexander F. Shchepetkin === ! ! ! This routine computes various horizontal metric terms. ! ! ! !======================================================================= ! implicit none PRIVATE PUBLIC :: metrics CONTAINS ! !*********************************************************************** SUBROUTINE metrics (ng, tile, model) !*********************************************************************** ! USE mod_param USE mod_grid USE mod_ocean USE mod_stepping ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model ! ! Local variable declarations. ! integer :: IminS, ImaxS, JminS, JmaxS integer :: LBi, UBi, LBj, UBj, LBij, UBij ! ! Set horizontal starting and ending indices for automatic private storage ! arrays. ! IminS=BOUNDS(ng)%Istr(tile)-3 ImaxS=BOUNDS(ng)%Iend(tile)+3 JminS=BOUNDS(ng)%Jstr(tile)-3 JmaxS=BOUNDS(ng)%Jend(tile)+3 ! ! Determine array lower and upper bounds in the I- and J-directions. ! LBi=BOUNDS(ng)%LBi(tile) UBi=BOUNDS(ng)%UBi(tile) LBj=BOUNDS(ng)%LBj(tile) UBj=BOUNDS(ng)%UBj(tile) ! ! Set array lower and upper bounds for MIN(I,J)- and MAX(I,J)-directions. ! LBij=BOUNDS(ng)%LBij UBij=BOUNDS(ng)%UBij ! CALL metrics_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp(ng), nnew(ng), & & GRID(ng) % f, & & GRID(ng) % h, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % pmask, & & GRID(ng) % rmask, & & GRID(ng) % angler, & & GRID(ng) % CosAngler, & & GRID(ng) % SinAngler, & & GRID(ng) % Hz, & & GRID(ng) % z_r, & & GRID(ng) % z_w, & & GRID(ng) % om_p, & & GRID(ng) % om_r, & & GRID(ng) % om_u, & & GRID(ng) % om_v, & & GRID(ng) % on_p, & & GRID(ng) % on_r, & & GRID(ng) % on_u, & & GRID(ng) % on_v, & & GRID(ng) % fomn, & & GRID(ng) % omn, & & GRID(ng) % pnom_p, & & GRID(ng) % pnom_r, & & GRID(ng) % pnom_u, & & GRID(ng) % pnom_v, & & GRID(ng) % pmon_p, & & GRID(ng) % pmon_r, & & GRID(ng) % pmon_u, & & GRID(ng) % pmon_v) RETURN END SUBROUTINE metrics ! !*********************************************************************** SUBROUTINE metrics_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, & & f, h, pm, pn, & & pmask, rmask, & & angler, CosAngler, SinAngler, & & Hz, z_r, z_w, & & om_p, om_r, om_u, om_v, & & on_p, on_r, on_u, on_v, & & fomn, omn, & & pnom_p, pnom_r, pnom_u, pnom_v, & & pmon_p, pmon_r, pmon_u, pmon_v) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_iounits USE mod_ncparam USE mod_scalars USE mod_iounits ! USE distribute_mod, ONLY : mp_reduce USE mp_exchange_mod, ONLY : mp_exchange2d USE set_depth_mod, ONLY : set_depth_tile implicit none ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nstp, nnew real(r8), intent(in) :: f(LBi:,LBj:) real(r8), intent(in) :: pm(LBi:,LBj:) real(r8), intent(in) :: pn(LBi:,LBj:) real(r8), intent(inout) :: h(LBi:,LBj:) real(r8), intent(inout) :: pmask(LBi:,LBj:) real(r8), intent(inout) :: rmask(LBi:,LBj:) real(r8), intent(inout) :: angler(LBi:,LBj:) real(r8), intent(out) :: om_p(LBi:,LBj:) real(r8), intent(out) :: om_r(LBi:,LBj:) real(r8), intent(out) :: om_u(LBi:,LBj:) real(r8), intent(out) :: om_v(LBi:,LBj:) real(r8), intent(out) :: on_p(LBi:,LBj:) real(r8), intent(out) :: on_r(LBi:,LBj:) real(r8), intent(out) :: on_u(LBi:,LBj:) real(r8), intent(out) :: on_v(LBi:,LBj:) real(r8), intent(out) :: fomn(LBi:,LBj:) real(r8), intent(out) :: omn(LBi:,LBj:) real(r8), intent(out) :: pnom_p(LBi:,LBj:) real(r8), intent(out) :: pnom_r(LBi:,LBj:) real(r8), intent(out) :: pnom_u(LBi:,LBj:) real(r8), intent(out) :: pnom_v(LBi:,LBj:) real(r8), intent(out) :: pmon_p(LBi:,LBj:) real(r8), intent(out) :: pmon_r(LBi:,LBj:) real(r8), intent(out) :: pmon_u(LBi:,LBj:) real(r8), intent(out) :: pmon_v(LBi:,LBj:) real(r8), intent(out) :: CosAngler(LBi:,LBj:) real(r8), intent(out) :: SinAngler(LBi:,LBj:) real(r8), intent(out) :: Hz(LBi:,LBj:,:) real(r8), intent(out) :: z_r(LBi:,LBj:,:) real(r8), intent(out) :: z_w(LBi:,LBj:,0:) ! ! Local variable declarations. ! logical :: EWperiodic=.FALSE. logical :: NSperiodic=.FALSE. integer :: NSUB, bry, i, is, j, k, rec real(r8), parameter :: Large = 1.0E+20_r8 real(r8) :: cff, cff1, cff2 real(r8) :: my_DXmax, my_DXmin, my_DYmax, my_DYmin real(r8) :: my_DZmax, my_DZmin real(r8) :: my_Cu_Cor, my_Cu_max, my_Cu_min, my_grdmax real(r8), dimension(14) :: buffer character (len=3), dimension(14) :: op_handle real(r8), dimension(LBi:UBi,LBj:UBj) :: A2d ! !----------------------------------------------------------------------- ! Set lower and upper tile bounds and staggered variables bounds for ! this horizontal domain partition. Notice that if tile=-1, it will ! set the values for the global grid. !----------------------------------------------------------------------- ! integer :: Istr, IstrR, IstrT, IstrU, Iend, IendR, IendT integer :: Jstr, JstrR, JstrT, JstrV, Jend, JendR, JendT ! Istr =BOUNDS(ng)%Istr (tile) IstrR=BOUNDS(ng)%IstrR(tile) IstrT=BOUNDS(ng)%IstrT(tile) IstrU=BOUNDS(ng)%IstrU(tile) Iend =BOUNDS(ng)%Iend (tile) IendR=BOUNDS(ng)%IendR(tile) IendT=BOUNDS(ng)%IendT(tile) Jstr =BOUNDS(ng)%Jstr (tile) JstrR=BOUNDS(ng)%JstrR(tile) JstrT=BOUNDS(ng)%JstrT(tile) JstrV=BOUNDS(ng)%JstrV(tile) Jend =BOUNDS(ng)%Jend (tile) JendR=BOUNDS(ng)%JendR(tile) JendT=BOUNDS(ng)%JendT(tile) ! !----------------------------------------------------------------------- ! Compute 1/m, 1/n, 1/mn, and f/mn at horizontal RHO-points. !----------------------------------------------------------------------- ! DO j=JstrR,JendR DO i=IstrR,IendR om_r(i,j)=1.0_r8/pm(i,j) on_r(i,j)=1.0_r8/pn(i,j) omn(i,j)=1.0_r8/(pm(i,j)*pn(i,j)) fomn(i,j)=f(i,j)*omn(i,j) END DO END DO CALL mp_exchange2d (ng, tile, model, 4, & & LBi, UBi, LBj, UBj, & & NghostPoints, EWperiodic, NSperiodic, & & om_r, on_r, omn, fomn) ! !----------------------------------------------------------------------- ! Compute n/m, and m/n at horizontal RHO-points. !----------------------------------------------------------------------- ! DO j=JstrR,JendR DO i=IstrR,IendR pnom_r(i,j)=pn(i,j)/pm(i,j) pmon_r(i,j)=pm(i,j)/pn(i,j) END DO END DO CALL mp_exchange2d (ng, tile, model, 2, & & LBi, UBi, LBj, UBj, & & NghostPoints, EWperiodic, NSperiodic, & & pnom_r, pmon_r) ! !----------------------------------------------------------------------- ! Compute m/n, 1/m, and 1/n at horizontal U-points. !----------------------------------------------------------------------- ! DO j=JstrR,JendR DO i=Istr,IendR pmon_u(i,j)=(pm(i-1,j)+pm(i,j))/(pn(i-1,j)+pn(i,j)) pnom_u(i,j)=(pn(i-1,j)+pn(i,j))/(pm(i-1,j)+pm(i,j)) om_u(i,j)=2.0_r8/(pm(i-1,j)+pm(i,j)) on_u(i,j)=2.0_r8/(pn(i-1,j)+pn(i,j)) END DO END DO CALL mp_exchange2d (ng, tile, model, 4, & & LBi, UBi, LBj, UBj, & & NghostPoints, EWperiodic, NSperiodic, & & pmon_u, pnom_u, om_u, on_u) ! !----------------------------------------------------------------------- ! Compute n/m, 1/m, and 1/m at horizontal V-points. !----------------------------------------------------------------------- ! DO j=Jstr,JendR DO i=IstrR,IendR pmon_v(i,j)=(pm(i,j-1)+pm(i,j))/(pn(i,j-1)+pn(i,j)) pnom_v(i,j)=(pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j)) om_v(i,j)=2.0_r8/(pm(i,j-1)+pm(i,j)) on_v(i,j)=2.0_r8/(pn(i,j-1)+pn(i,j)) END DO END DO CALL mp_exchange2d (ng, tile, model, 4, & & LBi, UBi, LBj, UBj, & & NghostPoints, EWperiodic, NSperiodic, & & pmon_v, pnom_v, om_v, on_v) ! !----------------------------------------------------------------------- ! Compute n/m and m/n at horizontal PSI-points. !----------------------------------------------------------------------- ! DO j=Jstr,JendR DO i=Istr,IendR pnom_p(i,j)=(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))/ & & (pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j)) pmon_p(i,j)=(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))/ & & (pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j)) om_p(i,j)=4.0_r8/(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j)) on_p(i,j)=4.0_r8/(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j)) END DO END DO CALL mp_exchange2d (ng, tile, model, 4, & & LBi, UBi, LBj, UBj, & & NghostPoints, EWperiodic, NSperiodic, & & pnom_p, pmon_p, om_p, on_p) ! !----------------------------------------------------------------------- ! Set slipperiness (no-slip) mask at PSI-points. !----------------------------------------------------------------------- ! ! Set no-slip boundary conditions on land-mask boundaries regardless of ! supplied value of gamma2. ! cff1=1.0_r8 ! computation of off-diagonal nonlinear terms cff2=2.0_r8 DO j=Jstr,Jend DO i=Istr,Iend IF ((rmask(i-1,j ).gt.0.5_r8).and. & & (rmask(i ,j ).gt.0.5_r8).and. & & (rmask(i-1,j-1).gt.0.5_r8).and. & & (rmask(i ,j-1).gt.0.5_r8)) THEN pmask(i,j)=1.0_r8 ELSE IF ((rmask(i-1,j ).lt.0.5_r8).and. & & (rmask(i ,j ).gt.0.5_r8).and. & & (rmask(i-1,j-1).gt.0.5_r8).and. & & (rmask(i ,j-1).gt.0.5_r8)) THEN pmask(i,j)=cff1 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. & & (rmask(i ,j ).lt.0.5_r8).and. & & (rmask(i-1,j-1).gt.0.5_r8).and. & & (rmask(i ,j-1).gt.0.5_r8)) THEN pmask(i,j)=cff1 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. & & (rmask(i ,j ).gt.0.5_r8).and. & & (rmask(i-1,j-1).lt.0.5_r8).and. & & (rmask(i ,j-1).gt.0.5_r8)) THEN pmask(i,j)=cff1 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. & & (rmask(i ,j ).gt.0.5_r8).and. & & (rmask(i-1,j-1).gt.0.5_r8).and. & & (rmask(i ,j-1).lt.0.5_r8)) THEN pmask(i,j)=cff1 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. & & (rmask(i ,j ).lt.0.5_r8).and. & & (rmask(i-1,j-1).gt.0.5_r8).and. & & (rmask(i ,j-1).lt.0.5_r8)) THEN pmask(i,j)=cff2 ELSE IF ((rmask(i-1,j ).lt.0.5_r8).and. & & (rmask(i ,j ).gt.0.5_r8).and. & & (rmask(i-1,j-1).lt.0.5_r8).and. & & (rmask(i ,j-1).gt.0.5_r8)) THEN pmask(i,j)=cff2 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. & & (rmask(i ,j ).gt.0.5_r8).and. & & (rmask(i-1,j-1).lt.0.5_r8).and. & & (rmask(i ,j-1).lt.0.5_r8)) THEN pmask(i,j)=cff2 ELSE IF ((rmask(i-1,j ).lt.0.5_r8).and. & & (rmask(i ,j ).lt.0.5_r8).and. & & (rmask(i-1,j-1).gt.0.5_r8).and. & & (rmask(i ,j-1).gt.0.5_r8)) THEN pmask(i,j)=cff2 ELSE pmask(i,j)=0.0_r8 END IF END DO END DO CALL mp_exchange2d (ng, tile, model, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, EWperiodic, NSperiodic, & & pmask) ! !----------------------------------------------------------------------- ! Compute cosine and sine of curvilinear rotation angle. !----------------------------------------------------------------------- ! DO j=JstrR,JendR DO i=IstrR,IendR CosAngler(i,j)=COS(angler(i,j)) SinAngler(i,j)=SIN(angler(i,j)) END DO END DO CALL mp_exchange2d (ng, tile, model, 2, & & LBi, UBi, LBj, UBj, & & NghostPoints, EWperiodic, NSperiodic, & & CosAngler, SinAngler) ! !----------------------------------------------------------------------- ! Compute minimum and maximum grid spacing. !----------------------------------------------------------------------- ! ! Compute time invariant depths (use zero free-surface). ! DO i=LBi,UBi DO j=LBj,UBj A2d(i,j)=0.0_r8 END DO END DO CALL set_depth_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, & & h, & & A2d, & & Hz, z_r, z_w) ! ! Compute grid spacing range. ! my_DXmin= Large my_DXmax=-Large my_DYmin= Large my_DYmax=-Large my_DZmin= Large my_DZmax=-Large my_grdmax=-Large DO j=JstrR,JendR DO i=IstrR,IendR cff=SQRT(om_r(i,j)*on_r(i,j)) my_DXmin=MIN(my_DXmin,om_r(i,j)) my_DXmax=MAX(my_DXmax,om_r(i,j)) my_DYmin=MIN(my_DYmin,on_r(i,j)) my_DYmax=MAX(my_DYmax,on_r(i,j)) my_grdmax=MAX(my_grdmax,cff) DO k=1,N(ng) my_DZmin=MIN(my_DZmin,Hz(i,j,k)) my_DZmax=MAX(my_DZmax,Hz(i,j,k)) END DO END DO END DO ! !----------------------------------------------------------------------- ! Compute Courant number. !----------------------------------------------------------------------- ! ! The Courant number is defined as: ! ! Cu = c * dt * SQRT (1/dx^2 + 1/dy^2) ! ! where c=SQRT(g*h) is phase speed for barotropic mode, and dx, dy ! are grid spacing in each direction. ! my_Cu_min= Large my_Cu_max=-Large my_Cu_Cor=-Large DO j=JstrR,JendR DO i=IstrR,IendR IF (rmask(i,j).gt.0.0_r8) THEN cff=dtfast(ng)* & & SQRT(g*ABS(h(i,j))*(pm(i,j)*pm(i,j)+pn(i,j)*pn(i,j))) my_Cu_min=MIN(my_Cu_min,cff) my_Cu_max=MAX(my_Cu_max,cff) cff=dt(ng)*ABS(f(i,j)) my_Cu_Cor=MAX(my_Cu_Cor,cff) END IF END DO END DO ! !----------------------------------------------------------------------- ! Perform global reductions. !----------------------------------------------------------------------- ! IF ((Istr.eq.1).and.(Jstr.eq.1).and. & & (Iend.eq.Lm(ng)).and.(Jend.eq.Mm(ng))) THEN NSUB=1 ! non-tiled application ELSE NSUB=NtileX(ng)*NtileE(ng) ! tiled application END IF IF (tile_count.eq.0) THEN Cu_min=my_Cu_min Cu_max=my_Cu_max Cu_Cor=my_Cu_Cor grdmax(ng)=my_grdmax DXmin(ng)=my_DXmin DXmax(ng)=my_DXmax DYmin(ng)=my_DYmin DYmax(ng)=my_DYmax DZmin(ng)=my_DZmin DZmax(ng)=my_DZmax ELSE Cu_min=MIN(Cu_min,my_Cu_min) Cu_max=MAX(Cu_max,my_Cu_max) Cu_Cor=MAX(Cu_Cor,my_Cu_Cor) grdmax(ng)=MAX(grdmax(ng),my_grdmax) DXmin(ng)=MIN(DXmin(ng),my_DXmin) DXmax(ng)=MAX(DXmax(ng),my_DXmax) DYmin(ng)=MIN(DYmin(ng),my_DYmin) DYmax(ng)=MAX(DYmax(ng),my_DYmax) DZmin(ng)=MIN(DZmin(ng),my_DZmin) DZmax(ng)=MAX(DZmax(ng),my_DZmax) END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 buffer(1)=Cu_min op_handle(1)='MIN' buffer(2)=Cu_max op_handle(2)='MAX' buffer(3)=Cu_Cor op_handle(3)='MAX' buffer(4)=grdmax(ng) op_handle(4)='MAX' buffer(5)=DXmin(ng) op_handle(5)='MIN' buffer(6)=DXmax(ng) op_handle(6)='MAX' buffer(7)=DYmin(ng) op_handle(7)='MIN' buffer(8)=DYmax(ng) op_handle(8)='MAX' buffer(9)=DZmin(ng) op_handle(9)='MIN' buffer(10)=DZmax(ng) op_handle(10)='MAX' buffer(11)=0.0_r8 op_handle(11)='MIN' buffer(12)=0.0_r8 op_handle(12)='MAX' buffer(13)=0.0_r8 op_handle(13)='MIN' buffer(14)=0.0_r8 op_handle(14)='MAX' CALL mp_reduce (ng, model, 14, buffer, op_handle) Cu_min=buffer(1) Cu_max=buffer(2) Cu_Cor=buffer(3) grdmax(ng)=buffer(4) DXmin(ng)=buffer(5) DXmax(ng)=buffer(6) DYmin(ng)=buffer(7) DYmax(ng)=buffer(8) DZmin(ng)=buffer(9) DZmax(ng)=buffer(10) IF (Master.and.LwrtInfo(ng)) THEN WRITE(stdout,10) DXmin(ng)/1000.0_r8, DXmax(ng)/1000.0_r8, & & DYmin(ng)/1000.0_r8, DYmax(ng)/1000.0_r8 10 FORMAT (/,' Minimum X-grid spacing, DXmin = ',1pe15.8,' km', & & /,' Maximum X-grid spacing, DXmax = ',1pe15.8,' km', & & /,' Minimum Y-grid spacing, DYmin = ',1pe15.8,' km', & & /,' Maximum Y-grid spacing, DYmax = ',1pe15.8,' km') WRITE(stdout,20) DZmin(ng), DZmax(ng) 20 FORMAT (' Minimum Z-grid spacing, DZmin = ',1pe15.8,' m',/, & & ' Maximum Z-grid spacing, DZmax = ',1pe15.8,' m') WRITE (stdout,30) Cu_min, Cu_max, Cu_Cor 30 FORMAT (/,' Minimum barotropic Courant Number = ', 1pe15.8,/, & & ' Maximum barotropic Courant Number = ', 1pe15.8,/, & & ' Maximum Coriolis Courant Number = ', 1pe15.8,/) END IF END IF RETURN END SUBROUTINE metrics_tile END MODULE metrics_mod