I have a setup with a spherical grid being forced by restoring buoyancy at the surface. This set up runs great.
For my future runs, because my focus is on the eastern boundary (EB), I want to set up the grid such that zonal resolution is high near EB and becomes coarse as we go towards the west. Not being familiar with mosaic, I decided to meddle around with MOM_grid_initialization.F90
to achieve this.
subroutine set_grid_metrics_spherical_nonuniformx(G, param_file)
type(ocean_grid_type), intent(inout) :: G
type(param_file_type), intent(in) :: param_file
! Arguments:
! (inout) G - The ocean's grid structure.
! (in) param_file - A structure indicating the open file to parse for
! model parameter values.
! Calculate the values of the metric terms that might be used
! and save them in arrays.
! Within this subroutine, the x- and y- grid spacings and their
! inverses and the cell areas centered on h, q, u, and v points are
! calculated, as are the geographic locations of each of these 4
! sets of points.
real :: PI, PI_180! PI = 3.1415926... as 4*atan(1)
integer :: i, j, isd, ied, jsd, jed
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
integer :: i_offset, j_offset, N
real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB)
real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB)
real :: dLon,dLat,latitude,longitude,dL_di, dLonmin, dLonmax
character(len=48) :: mod = "MOM_grid_init set_grid_metrics_spherical_nonuniformx"
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
i_offset = G%isd_global - isd; j_offset = G%jsd_global - jsd
call callTree_enter("set_grid_metrics_spherical_nonuniformx(), MOM_grid_initialize.F90")
! Calculate the values of the metric terms that might be used
! and save them in arrays.
PI = 4.0*atan(1.0); PI_180 = atan(1.0)/45.
call get_param(param_file, mod, "AXIS_UNITS", G%axis_units, default="degrees")
if (trim(G%axis_units) == "") G%axis_units = "degrees"
if (trim(G%axis_units) .ne. "degrees") call MOM_error(FATAL, &
"MOM_grid_init.F90, set_grid_metrics_simple_spherical_nonuniformx: "// &
"axis_units must be degrees")
call get_param(param_file, mod, "SOUTHLAT", G%south_lat, &
"The southern latitude of the domain.", units="degrees", &
fail_if_missing=.true.)
call get_param(param_file, mod, "LENLAT", G%len_lat, &
"The latitudinal length of the domain.", units="degrees", &
fail_if_missing=.true.)
call get_param(param_file, mod, "WESTLON", G%west_lon, &
"The western longitude of the domain.", units="degrees", &
default=0.0)
call get_param(param_file, mod, "LENLON", G%len_lon, &
"The longitudinal length of the domain.", units="degrees", &
fail_if_missing=.true.)
call get_param(param_file, mod, "RAD_EARTH", G%Rad_Earth, &
"The radius of the Earth.", units="m", default=6.378e6)
call get_param(param_file, mod, "DLON_MIN", dLonmin, &
"The desired minimum grid size in x direction.", units="degrees", &
fail_if_missing=.true.)
N = G%Domain%niglobal
! dLon = G%len_lon/G%Domain%niglobal
! dLonmin is the resolution near the eastern boundary and dLonmax is the resolution
! near the western boundary. dLonmax has to be calculated based upon domain length,
! number of grid points and dLonmin.
dLonmax = (PI*G%len_lon - 2.0*(N)*dLonmin)/((N)*(PI-2.0))
dLat = G%len_lat/G%Domain%njglobal
do j=G%JsgB,G%JegB
latitude = G%south_lat + dLat*(REAL(J-(G%jsg-1)))
G%gridLatB(J) = MIN(MAX(latitude,-90.),90.)
enddo
do j=G%jsg,G%jeg
latitude = G%south_lat + dLat*(REAL(j-G%jsg)+0.5)
G%gridLatT(j) = MIN(MAX(latitude,-90.),90.)
enddo
do i=G%IsgB,G%IegB
dlon = (dLonmin-dLonmax)*sin(i*PI/2.0/REAL(N)) + dLonmax
if (i==G%IsgB) then
G%gridLonB(I) = G%west_lon + dLon*(REAL(I-(G%isg-1)))
else
G%gridLonB(I) = G%gridLonB(I-1) + dLon
endif
enddo
do i=G%isg,G%ieg
! dlon = (dLonmin-dLonmax)*sin(i*PI/2.0/REAL(N)) + dLonmax
if (i==G%isg) then
G%gridLonT(i) = (G%west_lon + G%gridLonB(I))/2.0
else
G%gridLonT(i) = (G%gridLonB(I) + G%gridLonB(I-1))/2.0
endif
enddo
do J=JsdB,JedB
latitude = G%south_lat + dLat* REAL(J+J_offset-(G%jsg-1))
grid_LatB(J) = MIN(MAX(latitude,-90.),90.)
enddo
do j=jsd,jed
latitude = G%south_lat + dLat*(REAL(j+J_offset-G%jsg)+0.5)
grid_LatT(j) = MIN(MAX(latitude,-90.),90.)
enddo
do I=IsdB,IedB
dlon = (dLonmin-dLonmax)*sin(I*PI/2.0/REAL(N)) + dLonmax
if (i==IsdB) then
grid_LonB(I) = G%west_lon + dLon*REAL(I+I_offset-(G%isg-1))
else
grid_LonB(I) = grid_LonB(I-1) + dLon
endif
enddo
do i=isd,ied
! dlon = (dLonmin-dLonmax)*sin(i*PI/2.0/REAL(N)) + dLonmax
if (i==isd) then
grid_LonT(i) = (G%west_lon + grid_LonB(I))/2.0
else
grid_LonT(i) = (grid_LonB(i) + grid_LonB(i-1))/2.0
endif
enddo
! dL_di = (G%len_lon * 4.0*atan(1.0)) / (180.0 * G%Domain%niglobal)
do J=JsdB,JedB ; do I=IsdB,IedB
G%geoLonBu(I,J) = grid_lonB(I)
G%geoLatBu(I,J) = grid_latB(J)
dlon = (dLonmin-dLonmax)*sin(i*PI/2.0/REAL(N)) + dLonmax
! The following line is needed to reproduce the solution from
! set_grid_metrics_mercator when used to generate a simple spherical grid.
! G%dxBu(I,J) = G%Rad_Earth * COS( G%geoLatBu(I,J)*PI_180 ) * dL_di
G%dxBu(I,J) = G%Rad_Earth * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 )
G%dyBu(I,J) = G%Rad_Earth * dLat*PI_180
G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J)
enddo; enddo
do J=JsdB,JedB ; do i=isd,ied
G%geoLonCv(i,J) = grid_LonT(i)
G%geoLatCv(i,J) = grid_latB(J)
dlon = (dLonmin-dLonmax)*sin(i*PI/2.0/REAL(N)) + dLonmax
! The following line is needed to reproduce the solution from
! set_grid_metrics_mercator when used to generate a simple spherical grid.
! G%dxCv(i,J) = G%Rad_Earth * COS( G%geoLatCv(i,J)*PI_180 ) * dL_di
G%dxCv(i,J) = G%Rad_Earth * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 )
G%dyCv(i,J) = G%Rad_Earth * dLat*PI_180
enddo; enddo
do j=jsd,jed ; do I=IsdB,IedB
G%geoLonCu(I,j) = grid_lonB(I)
G%geoLatCu(I,j) = grid_LatT(j)
dlon = (dLonmin-dLonmax)*sin(i*PI/2.0/REAL(N)) + dLonmax
! The following line is needed to reproduce the solution from
! set_grid_metrics_mercator when used to generate a simple spherical grid.
! G%dxCu(I,j) = G%Rad_Earth * COS( G%geoLatCu(I,j)*PI_180 ) * dL_di
G%dxCu(I,j) = G%Rad_Earth * dLon*PI_180 * COS( G%geoLatCu(i,J)*PI_180 )
G%dyCu(I,j) = G%Rad_Earth * dLat*PI_180
enddo; enddo
do j=jsd,jed ; do i=isd,ied
G%geoLonT(i,j) = grid_LonT(i)
G%geoLatT(i,j) = grid_LatT(j)
dlon = (dLonmin-dLonmax)*sin(i*PI/2.0/REAL(N)) + dLonmax
! The following line is needed to reproduce the solution from
! set_grid_metrics_mercator when used to generate a simple spherical grid.
! G%dxT(i,j) = G%Rad_Earth * COS( G%geoLatT(i,j)*PI_180 ) * dL_di
G%dxT(i,j) = G%Rad_Earth * dLon*PI_180 * COS( G%geoLatT(i,j)*PI_180 )
G%dyT(i,j) = G%Rad_Earth * dLat*PI_180
! latitude = G%geoLatCv(i,J)*PI_180 ! In radians
! dL_di = G%geoLatCv(i,max(jsd,J-1))*PI_180 ! In radians
! G%areaT(i,j) = Rad_Earth**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di))
G%areaT(i,j) = G%dxT(i,j) * G%dyT(i,j)
enddo; enddo
call callTree_leave("set_grid_metrics_spherical_nonuniformx()")
end subroutine set_grid_metrics_spherical_nonuniformx