mosaic/gradient.F90

   1 !***********************************************************************
   2 !*                   GNU Lesser General Public License
   3 !*
   4 !* This file is part of the GFDL Flexible Modeling System (FMS).
   5 !*
   6 !* FMS is free software: you can redistribute it and/or modify it under
   7 !* the terms of the GNU Lesser General Public License as published by
   8 !* the Free Software Foundation, either version 3 of the License, or (at
   9 !* your option) any later version.
  10 !*
  11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
  12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 !* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 !* for more details.
  15 !*
  16 !* You should have received a copy of the GNU Lesser General Public
  17 !* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18 !***********************************************************************
  19 !> @defgroup gradient_mod gradient_mod
  20 !> @ingroup mosaic
  21 !> @brief Implements some utility routines to calculate gradient.
  22 !> @author Zhi Liang
  23 !!
  24 !! Currently only gradient on cubic grid is implemented. Also a public interface
  25 !! is provided to calculate grid information needed to calculate gradient.
  26
  27 !> @addtogroup gradient_mod
  28 !> @{
  29 module gradient_mod
  30
  31 use mpp_mod,       only : mpp_error, FATAL
  32 use constants_mod, only : RADIUS
  33 use platform_mod
  34
  35 implicit none
  36 private
  37
  38
  39 public :: gradient_cubic
  40 public :: calc_cubic_grid_info
  41
  42 ! Include variable "version" to be written to log file.
  43 #include<file_version.h>
  44
  45 contains
  46
  47
  48 !#####################################################################
  49 !> Pin has halo size = 1.
  50 !! @param pin the size of pin    will be (nx+2,ny+2), T-cell center, with halo = 1
  51 !! @param the size of dx     will be (nx, ny+1),  N-cell center
  52 !! @param the size of dy     will be (nx+1, ny),  E-cell center
  53 !! @param the size of area   will be (nx, ny),    T-cell center.
  54 !! @param The size of edge_w will be (ny+1),      C-cell center
  55 !! @param The size of edge_e will be (ny+1),      C-cell center
  56 !! @param The size of edge_s will be (nx+1),      C-cell center
  57 !! @param The size of edge_n will be (nx+1),      C-cell center
  58 !! @param The size of en_n   will be (3,nx,ny+1), N-cell center
  59 !! @param The size of en_e   will be (3,nx+1,ny), E-cell center
  60 !! @param The size of vlon   will be (3,nx, ny)   T-cell center
  61 !! @param The size of vlat   will be (3,nx, ny),  T-cell center
  62 subroutine gradient_cubic(pin, dx, dy, area, edge_w, edge_e, edge_s, edge_n,    &
  63                           en_n, en_e, vlon, vlat, grad_x, grad_y, on_west_edge, &
  64                           on_east_edge, on_south_edge, on_north_edge)
  65
  66   real(r8_kind),    dimension(:,:  ), intent(in ) :: pin, dx, dy, area
  67   real(r8_kind),    dimension(:    ), intent(in ) :: edge_w, edge_e, edge_s, edge_n
  68   real(r8_kind),    dimension(:,:,:), intent(in ) :: en_n, en_e
  69   real(r8_kind),    dimension(:,:,:), intent(in ) :: vlon, vlat
  70   real(r8_kind),    dimension(:,:  ), intent(out) :: grad_x, grad_y
  71   logical,                   intent(in ) :: on_west_edge, on_east_edge, on_south_edge, on_north_edge
  72   integer :: nx, ny
  73
  74
  75   nx = size(grad_x,1)
  76   ny = size(grad_x,2)
  77
  78   if(size(pin,1) .NE. nx+2 .OR. size(pin,2) .NE. ny+2)call mpp_error(FATAL, &
  79      &  "gradient_mod:size of pin should be (nx+2, ny+2)")
  80   if(size(dx,1) .NE. nx .OR. size(dx,2) .NE. ny+1 ) call mpp_error(FATAL, &
  81      &  "gradient_mod: size of dx should be (nx,ny+1)")
  82   if(size(dy,1) .NE. nx+1 .OR. size(dy,2) .NE. ny ) call mpp_error(FATAL, &
  83      &  "gradient_mod: size of dy should be (nx+1,ny)")
  84   if(size(area,1) .NE. nx .OR. size(area,2) .NE. ny ) call mpp_error(FATAL, &
  85      &  "gradient_mod: size of area should be (nx,ny)")
  86   if(size(vlon,1) .NE. 3 .OR. size(vlon,2) .NE. nx .OR. size(vlon,3) .NE. ny) &
  87           call mpp_error(FATAL, "gradient_mod: size of vlon should be (3,nx,ny)")
  88   if(size(vlat,1) .NE. 3 .OR. size(vlat,2) .NE. nx .OR. size(vlat,3) .NE. ny) &
  89           call mpp_error(FATAL, "gradient_mod: size of vlat should be (3,nx,ny)")
  90   if(size(edge_w) .NE. ny+1) call mpp_error(FATAL, "gradient_mod: size of edge_w should be (ny+1)")
  91   if(size(edge_e) .NE. ny+1) call mpp_error(FATAL, "gradient_mod: size of edge_e should be (ny+1)")
  92   if(size(edge_s) .NE. nx+1) call mpp_error(FATAL, "gradient_mod: size of edge_s should be (nx+1)")
  93   if(size(edge_n) .NE. nx+1) call mpp_error(FATAL, "gradient_mod: size of edge_n should be (nx+1)")
  94   if(size(en_n,1) .NE. 3 .OR. size(en_n,2) .NE. nx .OR.  size(en_n,3) .NE. ny+1 ) &
  95        call mpp_error(FATAL, "gradient_mod:size of en_n should be (3, nx, ny+1)")
  96   if(size(en_e,1) .NE. 3 .OR. size(en_e,2) .NE. nx+1 .OR.  size(en_e,3) .NE. ny ) &
  97        call mpp_error(FATAL, "gradient_mod:size of en_e should be (3, nx+1, ny)")
  98
  99   call grad_c2l(nx, ny, pin, dx, dy, area, edge_w, edge_e, edge_s, edge_n, en_n, en_e, vlon, vlat, &
 100                 grad_x, grad_y, on_west_edge, on_east_edge, on_south_edge, on_north_edge)
 101
 102   return
 103
 104 end subroutine gradient_cubic
 105
 106
 107 subroutine calc_cubic_grid_info(xt, yt, xc, yc, dx, dy, area, edge_w, edge_e, edge_s, edge_n, &
 108                            en_n, en_e, vlon, vlat, on_west_edge, on_east_edge, on_south_edge, on_north_edge )
 109   real(r8_kind),    dimension(:,:  ), intent(in ) :: xt, yt, xc, yc
 110   real(r8_kind),    dimension(:,:  ), intent(out) :: dx, dy, area
 111   real(r8_kind),    dimension(:    ), intent(out) :: edge_w, edge_e, edge_s, edge_n
 112   real(r8_kind),    dimension(:,:,:), intent(out) :: en_n, en_e
 113   real(r8_kind),    dimension(:,:,:), intent(out) :: vlon, vlat
 114   logical,                   intent(in ) :: on_west_edge, on_east_edge, on_south_edge, on_north_edge
 115   integer :: nx, ny, nxp, nyp
 116
 117
 118   nx  = size(area,1)
 119   ny  = size(area,2)
 120   nxp = nx+1
 121   nyp = ny+1
 122
 123   if(size(xt,1) .NE. nx+2 .OR. size(xt,2) .NE. ny+2 ) call mpp_error(FATAL, &
 124      &  "gradient_mod: size of xt should be (nx+2,ny+2)")
 125   if(size(yt,1) .NE. nx+2 .OR. size(yt,2) .NE. ny+2 ) call mpp_error(FATAL, &
 126      &  "gradient_mod: size of yt should be (nx+2,ny+2)")
 127   if(size(xc,1) .NE. nxp .OR. size(xc,2) .NE. nyp ) call mpp_error(FATAL, &
 128      &  "gradient_mod: size of xc should be (nx+1,ny+1)")
 129   if(size(yc,1) .NE. nxp .OR. size(yc,2) .NE. nyp ) call mpp_error(FATAL, &
 130      &  "gradient_mod: size of yc should be (nx+1,ny+1)")
 131   if(size(dx,1) .NE. nx .OR. size(dx,2) .NE. nyp ) call mpp_error(FATAL, &
 132      &  "gradient_mod: size of dx should be (nx,ny+1)")
 133   if(size(dy,1) .NE. nxp .OR. size(dy,2) .NE. ny ) call mpp_error(FATAL, &
 134      &  "gradient_mod: size of dy should be (nx+1,ny)")
 135   if(size(area,1) .NE. nx .OR. size(area,2) .NE. ny ) call mpp_error(FATAL, &
 136      &  "gradient_mod: size of area should be (nx,ny)")
 137   if(size(vlon,1) .NE. 3 .OR. size(vlon,2) .NE. nx .OR. size(vlon,3) .NE. ny) &
 138           call mpp_error(FATAL, "gradient_mod: size of vlon should be (3,nx,ny)")
 139   if(size(vlat,1) .NE. 3 .OR. size(vlat,2) .NE. nx .OR. size(vlat,3) .NE. ny) &
 140           call mpp_error(FATAL, "gradient_mod: size of vlat should be (3,nx,ny)")
 141   if(size(edge_w) .NE. ny+1) call mpp_error(FATAL, "gradient_mod: size of edge_w should be (ny-1)")
 142   if(size(edge_e) .NE. ny+1) call mpp_error(FATAL, "gradient_mod: size of edge_e should be (ny-1)")
 143   if(size(edge_s) .NE. nx+1) call mpp_error(FATAL, "gradient_mod: size of edge_s should be (nx-1)")
 144   if(size(edge_n) .NE. nx+1) call mpp_error(FATAL, "gradient_mod: size of edge_n should be (nx-1)")
 145   if(size(en_n,1) .NE. 3 .OR. size(en_n,2) .NE. nx .OR.  size(en_n,3) .NE. nyp ) &
 146        call mpp_error(FATAL, "gradient_mod:size of en_n should be (3, nx, ny+1)")
 147   if(size(en_e,1) .NE. 3 .OR. size(en_e,2) .NE. nxp .OR.  size(en_e,3) .NE. ny ) &
 148        call mpp_error(FATAL, "gradient_mod:size of en_e should be (3, nx+1, ny)")
 149
 150
 151   call calc_c2l_grid_info(nx, ny, xt, yt, xc, yc, dx, dy, area, edge_w, edge_e, edge_s, edge_n, &
 152                           en_n, en_e, vlon, vlat, on_west_edge, on_east_edge, on_south_edge, on_north_edge )
 153
 154
 155   return
 156
 157 end subroutine calc_cubic_grid_info
 158
 159 end module gradient_mod
 160 !> @}
 161 ! close documentation grouping