Line data Source code
1 : module edyn_geogrid
2 : !
3 : ! Global geographic grid.
4 : ! See sub set_geogrid (edyn_init.F90)
5 : !
6 : use shr_kind_mod, only: r8 => shr_kind_r8 ! 8-byte reals
7 : use cam_logfile, only: iulog
8 : use cam_abortutils, only: endrun
9 :
10 : implicit none
11 : private
12 : save
13 :
14 : integer, public, protected :: & ! dimensions
15 : nlat, & ! number of latitudes
16 : nlon, & ! number of longitudes
17 : nlev, & ! number of midpoint levels
18 : nilev, & ! number of interface levels
19 : npes ! number of PEs in geogrid
20 :
21 : real(r8), public, protected, allocatable , dimension(:) :: & ! coordinate vars
22 : glat, & ! latitude coordinates (degrees)
23 : glon, & ! longitude coordinates (degrees)
24 : ylatg, & ! latitudes (radians)
25 : ylong, & ! longitudes (radians)
26 : zlev, & ! midpoint vertical coordinates
27 : zilev ! interface vertical coordinates
28 :
29 : real(r8), public, allocatable, protected :: cs(:) ! cos(phi) (0:nlat+1)
30 :
31 : integer, public, protected :: & ! model independent (set by sub get_geogrid)
32 : nlonp1, & ! nlon+1
33 : nlonp2, & ! nlon+2
34 : nlatp1 ! nlat+1
35 :
36 : real(r8), public, protected :: dphi
37 : real(r8), public, protected :: dlamda
38 : !
39 : ! Using p0 in microbars, as in TIEGCM.
40 : real(r8), parameter, public :: p0 = 5.0e-4_r8 ! standard pressure (microbars)
41 :
42 : integer, public, protected :: & ! model dependent (set by subs read_tgcm, read_waccm)
43 : jspole, & ! latitude index to geographic south pole
44 : jnpole ! latitude index to geographic north pole
45 :
46 : ! set_geogrid sets up a distributed finite-volume lat/lon grid
47 : public :: set_geogrid
48 :
49 : logical :: debug = .false. ! set true for prints to stdout at each call
50 :
51 : contains
52 :
53 : !-----------------------------------------------------------------------
54 768 : subroutine set_geogrid(nlon_g, nlat_g, nlev_in, npes_in, iam, pres_mid_in, pres_edge_in, min_lat_pe_in)
55 : use shr_const_mod, only: pi => shr_const_pi
56 : use edyn_params, only: kbotdyn, pbotdyn
57 : use edyn_mpi, only: mp_distribute_geo
58 : use spmd_utils, only: masterproc
59 : use edyn_maggrid, only: nmlat
60 :
61 : ! Dummy Args
62 : integer, intent(in) :: nlon_g ! Global num longitudes
63 : integer, intent(in) :: nlat_g ! Global num latitudes
64 : integer, intent(in) :: nlev_in ! Num levels
65 : integer, intent(in) :: npes_in
66 : integer, intent(in) :: iam
67 : real(r8), intent(in) :: pres_mid_in(:)
68 : real(r8), intent(in) :: pres_edge_in(:)
69 : integer, optional, intent(in) :: min_lat_pe_in ! Min # lats / PE
70 : !
71 : ! Local:
72 : integer :: latind, lonind, js, k
73 : integer :: lon_beg, lon_end, lat_beg, lat_end
74 : integer :: lons_per_task, lats_per_task
75 : integer :: lons_overflow, lats_overflow
76 : integer :: ntasks_lat, ntasks_lon
77 : integer :: task_cnt, i,j
78 : integer :: minlats_per_pe
79 : integer :: ierr
80 : real(r8) :: phi
81 : real(r8) :: delta ! Coordinate spacing
82 : real(r8), parameter :: eps = 1.e-6_r8
83 :
84 1536 : real(r8) :: pmid(nlev_in)
85 :
86 768 : nlon = nlon_g
87 768 : nlat = nlat_g
88 768 : nlev = nlev_in
89 768 : npes = npes_in
90 :
91 768 : nilev = nlev+1
92 :
93 768 : nlonp1 = nlon + 1
94 768 : nlonp2 = nlon + 2
95 768 : nlatp1 = nlat + 1
96 :
97 768 : jspole = 1
98 768 : jnpole = nlat
99 :
100 768 : if (present(min_lat_pe_in)) then
101 0 : minlats_per_pe = min_lat_pe_in
102 : else
103 : minlats_per_pe = 2
104 : end if
105 :
106 768 : dphi = pi / real(nlat,r8)
107 768 : dlamda = 2._r8*pi / real(nlon,r8)
108 :
109 : !
110 : ! Allocate coordinate variables:
111 : !
112 2304 : allocate(glon(nlon))
113 2304 : allocate(glat(nlat))
114 : !
115 : ! Create a finite-volume coordinate grid (in degrees)
116 : !
117 768 : delta = 360.0_r8 / real(nlon, r8)
118 111360 : do lonind = 1, nlon
119 111360 : glon(lonind) = -180.0_r8 + ((lonind - 1) * delta)
120 : end do
121 768 : delta = 180.0_r8 / real((nlat - 1), r8)
122 : ! Set the poles exactly (they might be checked later)
123 768 : glat(1) = -90.0_r8
124 768 : glat(nlat) = 90.0_r8
125 72960 : do latind = 2, nlat - 1
126 72960 : glat(latind) = -90.0_r8 + ((latind - 1) * delta)
127 : end do
128 :
129 768 : if (masterproc.and.debug) then
130 0 : write(iulog,*) 'set_geogrid glon : ',glon(:)
131 0 : write(iulog,*) 'set_geogrid glat : ',glat(:)
132 : end if
133 :
134 2304 : allocate(zlev(nlev))
135 2304 : allocate(zilev(nilev))
136 : !
137 : ! Hybrid-sigma levels from ref_pres module:
138 : !
139 100608 : zlev(:nlev) = pres_mid_in(:) ! midpoints vertical coord (top down)
140 101376 : zilev(:nilev) = pres_edge_in(:nilev) ! interfaces vertical coord
141 :
142 : ! do bottom up search for kbotdyn
143 100608 : pmid(:nlev) = zlev(nlev:1:-1)
144 47616 : kloop: do k = 1, nlev
145 47616 : if ( pmid(k) <= pbotdyn) then
146 768 : kbotdyn = k
147 768 : exit kloop
148 : end if
149 : end do kloop
150 768 : if ( kbotdyn < 1 ) then
151 0 : call endrun('set_geogrid: kbotdyn is not set')
152 : endif
153 768 : if (debug) then
154 0 : write(iulog,"('set_geogrid: kbotdyn=',i4,' pmid(kbotdyn)=',es12.4)") kbotdyn,pmid(kbotdyn)
155 : endif
156 :
157 : !
158 : ! Setup a decomposition for the geogrid
159 : !
160 : ! First, try using a 1-D latitude decomposition
161 :
162 9216 : do ntasks_lon = 1,nlon_g
163 9216 : ntasks_lat = npes/ntasks_lon
164 9216 : if ( minlats_per_pe*ntasks_lat<nmlat .and. minlats_per_pe*ntasks_lat<nlat_g .and. ntasks_lat*ntasks_lon==npes ) then
165 : exit
166 : endif
167 : end do
168 768 : if (masterproc) then
169 2 : write(iulog,'(a,3i6)') 'set_geogrid: npes,nlon_g,nlat_g: ',npes,nlon_g,nlat_g
170 2 : write(iulog,'(a,2i6)') 'set_geogrid: ntasks_lon,ntasks_lat (oplus and edyn grids): ',ntasks_lon,ntasks_lat
171 : endif
172 :
173 768 : if (ntasks_lat*ntasks_lon/=npes) then
174 0 : call endrun('set_geogrid: ntasks_lat*ntasks_lon/=npes')
175 : endif
176 :
177 : ! Now, figure the starting and ending coordinates
178 768 : lons_per_task = nlon / ntasks_lon
179 768 : lons_overflow = MOD(nlon, ntasks_lon)
180 768 : lats_per_task = nlat / ntasks_lat
181 768 : lats_overflow = MOD(nlat, ntasks_lat)
182 768 : lon_beg = 1
183 768 : lon_end = 0
184 768 : lat_beg = 1
185 768 : lat_end = 0
186 768 : task_cnt= 0
187 768 : if (iam<npes) then
188 12672 : jloop: do j = 0,ntasks_lat-1
189 12672 : lat_beg = lat_end + 1
190 12672 : lat_end = lat_beg + lats_per_task - 1
191 12672 : if (j<lats_overflow) then
192 0 : lat_end = lat_end + 1
193 : end if
194 12672 : lon_end = 0
195 160512 : do i = 0,ntasks_lon-1
196 147840 : lon_beg = lon_end + 1
197 147840 : lon_end = lon_beg + lons_per_task - 1
198 147840 : if (i<lons_overflow) then
199 0 : lon_end = lon_end + 1
200 : end if
201 147840 : task_cnt = task_cnt+1
202 159744 : if (task_cnt>iam) exit jloop
203 : end do
204 : enddo jloop
205 : endif
206 :
207 768 : call mp_distribute_geo(lon_beg, lon_end, lat_beg, lat_end, 1, nlev, ntasks_lon, ntasks_lat)
208 :
209 : !
210 : ! Set horizontal geographic grid in radians (for apex code):
211 : !
212 2304 : allocate(ylatg(nlat)) ! waccm grid includes poles
213 2304 : allocate(ylong(nlonp1)) ! single periodic point
214 768 : ylatg(1) = -pi/2._r8+eps ! south pole
215 768 : ylatg(nlat) = pi/2._r8-eps ! north pole
216 72960 : do latind = 2, nlat-1
217 72960 : ylatg(latind) = -0.5_r8*(pi-dphi)+real(latind-1,r8)*dphi
218 : end do
219 112128 : do lonind = 1, nlonp1
220 112128 : ylong(lonind) = -pi+real(lonind-1,r8)*dlamda
221 : end do
222 : !
223 : ! Calculate cosine of latitude
224 : !
225 2304 : allocate(cs(0:nlat+1))
226 768 : js = -(nlat/2)
227 74496 : do latind = 1, nlat
228 73728 : phi = (latind + js - .5_r8) * dphi
229 74496 : cs(latind) = cos(phi)
230 : end do
231 768 : cs(0) = -cs(1)
232 768 : cs(nlat+1) = -cs(nlat)
233 :
234 768 : end subroutine set_geogrid
235 :
236 : end module edyn_geogrid
|
