Line data Source code
1 : module oplus
2 : !
3 : ! Horizontally transport the O+ ion, adapted for WACCM-X from TIEGCM.
4 : ! Input O+ is received from WACCM physics/chemistry, transported O+
5 : ! (op_out and opnm_out) are passed back to chemistry.
6 : !
7 : ! B. Foster (foster@ucar.edu), May, 2015.
8 : !
9 : use shr_kind_mod, only: r8 => shr_kind_r8
10 : use shr_const_mod, only: shr_const_g ! gravitational constant (m/s^2)
11 : use cam_abortutils, only: endrun
12 : use cam_logfile, only: iulog
13 : use spmd_utils, only: masterproc
14 : use savefield_waccm, only: savefld_waccm ! save field to waccm history
15 : use edyn_geogrid, only: dphi, dlamda, cs, p0
16 : use getapex, only: bx, by, bz, bmod2 ! (0:nlonp1,jspole-1:jnpole+1)
17 : use edyn_params, only: Rearth ! Radius of Earth (cm)
18 : use time_manager, only: get_step_size, is_first_step, is_first_restart_step
19 : use edyn_mpi, only: array_ptr_type
20 : use infnan, only: nan, assignment(=)
21 :
22 : implicit none
23 : private
24 :
25 : public :: oplus_xport, oplus_init
26 : public :: kbot
27 :
28 : real(r8) :: pi, rtd
29 : !
30 : ! Constants in CGS:
31 : !
32 : real(r8),parameter :: boltz = 1.38E-16_r8 ! boltzman's constant (erg/kelvin)
33 : real(r8),parameter :: gask = 8.314e7_r8 ! gas constant (erg/mol)
34 : !
35 : ! Collision factor (tuneable) (see also local colfac in iondrag.F90)
36 : ! FIX: Collision factor colfac is set locally in iondrag.F90 and here.
37 : ! It should be in one location and shared between ionosphere and
38 : ! dpie_coupling.
39 : !
40 : real(r8),parameter :: colfac = 1.5_r8 ! see also iondrag.F90
41 : !
42 : ! Reciprocal of molecular mass (multiply is cheaper than divide)
43 : real(r8),parameter :: rmassinv_o2=1._r8/32._r8, rmassinv_o1=1._r8/16._r8, &
44 : rmassinv_n2=1._r8/28._r8
45 : real(r8),parameter :: rmass_op=16._r8
46 :
47 : real(r8) :: dzp ! delta zp (typically 0.5 from kbot to top)
48 : real(r8) :: grav_cm ! gravitational constant (cm/s^2)
49 : integer, protected :: kbot = -999 ! k-index corresponding to ~pbot
50 : real(r8),parameter :: pbot = 0.004_r8 ! Pa -- bottom of O+ transport (near 120 km)
51 : !
52 : ! The shapiro constant .03 is used for spatial smoothing of oplus,
53 : ! (shapiro is tuneable, and maybe should be a function of timestep size).
54 : ! dtsmooth and dtsmooth_div2 are used in the time smoothing.
55 : ! To turn off all smoothing here, set shapiro=0. and dtsmooth = 1.
56 : !
57 :
58 : real(r8),parameter :: &
59 : dtsmooth = 0.95_r8, & ! for time smoother
60 : dtsmooth_div2 = 0.5_r8*(1._r8-dtsmooth)
61 :
62 : real(r8) :: adiff_limiter
63 : real(r8) :: shapiro_const
64 : logical :: enforce_floor
65 : logical :: ring_polar_filter = .false.
66 : logical, parameter :: debug = .false.
67 :
68 : real(r8), allocatable :: expz(:) ! exp(-zp)
69 : real(r8), allocatable :: zp(:) ! log pressure (as in tiegcm lev(nlev))
70 :
71 : logical :: debug_hist
72 :
73 : contains
74 :
75 : !-----------------------------------------------------------------------
76 768 : subroutine oplus_init( adiff_limiter_in, shapiro_const_in, enforce_floor_in, ring_polar_filter_in, ionos_debug_hist )
77 :
78 : use cam_history, only : addfld, horiz_only
79 : use filter_module,only : filter_init
80 : use edyn_geogrid, only : nlev
81 :
82 : real(r8), intent(in) :: adiff_limiter_in
83 : real(r8), intent(in) :: shapiro_const_in
84 : logical , intent(in) :: enforce_floor_in
85 : logical , intent(in) :: ring_polar_filter_in
86 : logical , intent(in) :: ionos_debug_hist
87 :
88 768 : debug_hist = ionos_debug_hist
89 :
90 768 : shapiro_const = shapiro_const_in
91 768 : enforce_floor = enforce_floor_in
92 768 : adiff_limiter = adiff_limiter_in
93 768 : ring_polar_filter = ring_polar_filter_in
94 :
95 768 : call filter_init( ring_polar_filter )
96 :
97 1536 : call addfld ('op_dt' , (/ 'lev' /), 'I', ' ' ,'op_dt' , gridname='geo_grid')
98 1536 : call addfld ('amb_diff' , (/ 'lev' /), 'I', ' ' ,'amb_diff' , gridname='geo_grid')
99 1536 : call addfld ('dfield' , (/ 'lev' /), 'I', ' ' ,'dfield' , gridname='geo_grid')
100 1536 : call addfld ('dwind' , (/ 'lev' /), 'I', ' ' ,'dwind' , gridname='geo_grid')
101 :
102 768 : if (debug_hist) then
103 : !
104 : ! Save fields from oplus module:
105 : !
106 0 : call addfld ('OPLUS_Z' ,(/ 'lev' /), 'I', 'cm ','OPLUS_Z' , gridname='geo_grid')
107 0 : call addfld ('OPLUS_TN' ,(/ 'lev' /), 'I', 'deg K','OPLUS_TN' , gridname='geo_grid')
108 0 : call addfld ('OPLUS_TE' ,(/ 'lev' /), 'I', 'deg K','OPLUS_TE' , gridname='geo_grid')
109 0 : call addfld ('OPLUS_TI' ,(/ 'lev' /), 'I', 'deg K','OPLUS_TI' , gridname='geo_grid')
110 0 : call addfld ('OPLUS_UN' ,(/ 'lev' /), 'I', 'cm/s' ,'OPLUS_UN' , gridname='geo_grid')
111 0 : call addfld ('OPLUS_VN' ,(/ 'lev' /), 'I', 'cm/s' ,'OPLUS_VN' , gridname='geo_grid')
112 0 : call addfld ('OPLUS_OM' ,(/ 'lev' /), 'I', 'Pa/s' ,'OPLUS_OM' , gridname='geo_grid')
113 0 : call addfld ('OPLUS_O2' ,(/ 'lev' /), 'I', 'mmr' ,'OPLUS_O2' , gridname='geo_grid')
114 0 : call addfld ('OPLUS_O1' ,(/ 'lev' /), 'I', 'mmr' ,'OPLUS_O1' , gridname='geo_grid')
115 :
116 0 : call addfld ('OPLUS_N2' ,(/ 'lev' /), 'I', 'mmr' ,'OPLUS_N2' , gridname='geo_grid')
117 0 : call addfld ('OPLUS_OP' ,(/ 'lev' /), 'I', 'cm^3' ,'OPLUS_OP' , gridname='geo_grid')
118 0 : call addfld ('OPLUS_UI' ,(/ 'lev' /), 'I', 'm/s' ,'OPLUS_UI' , gridname='geo_grid')
119 0 : call addfld ('OPLUS_VI' ,(/ 'lev' /), 'I', 'm/s' ,'OPLUS_VI' , gridname='geo_grid')
120 0 : call addfld ('OPLUS_WI' ,(/ 'lev' /), 'I', 'm/s' ,'OPLUS_WI' , gridname='geo_grid')
121 0 : call addfld ('OPLUS_MBAR' ,(/ 'lev' /), 'I', ' ' ,'OPLUS_MBAR' , gridname='geo_grid')
122 0 : call addfld ('OPLUS_TR' ,(/ 'lev' /), 'I', ' ' ,'OPLUS_TR' , gridname='geo_grid')
123 0 : call addfld ('OPLUS_TP0' ,(/ 'lev' /), 'I', ' ' ,'OPLUS_TP0' , gridname='geo_grid')
124 0 : call addfld ('OPLUS_TP1' ,(/ 'lev' /), 'I', ' ' ,'OPLUS_TP1' , gridname='geo_grid')
125 0 : call addfld ('OPLUS_DJ' ,(/ 'lev' /), 'I', ' ' ,'OPLUS_DJ' , gridname='geo_grid')
126 0 : call addfld ('OPLUS_HJ' ,(/ 'lev' /), 'I', ' ' ,'OPLUS_HJ' , gridname='geo_grid')
127 0 : call addfld ('OPLUS_BVEL' ,(/ 'lev' /), 'I', ' ' ,'OPLUS_BVEL' , gridname='geo_grid')
128 0 : call addfld ('OPLUS_DIFFJ',(/ 'lev' /), 'I', ' ' ,'OPLUS_DIFFJ' , gridname='geo_grid')
129 0 : call addfld ('OPLUS_OPNM' ,(/ 'lev' /), 'I', ' ' ,'OPLUS_OPNM' , gridname='geo_grid')
130 0 : call addfld ('OPNM_SMOOTH',(/ 'lev' /), 'I', ' ' ,'OPNM_SMOOTH' , gridname='geo_grid')
131 0 : call addfld ('BDOTDH_OP' ,(/ 'lev' /), 'I', ' ' ,'BDOTDH_OP' , gridname='geo_grid')
132 0 : call addfld ('BDOTDH_OPJ' ,(/ 'lev' /), 'I', ' ' ,'BDOTDH_OPJ' , gridname='geo_grid')
133 0 : call addfld ('BDOTDH_DIFF',(/ 'lev' /), 'I', ' ' ,'BDOTDH_DIFF' , gridname='geo_grid')
134 0 : call addfld ('BDZDVB_OP' ,(/ 'lev' /), 'I', ' ' ,'BDZDVB_OP' , gridname='geo_grid')
135 0 : call addfld ('EXPLICIT0' ,(/ 'lev' /), 'I', ' ' ,'EXPLICIT0' , gridname='geo_grid')
136 :
137 0 : call addfld ('EXPLICITa' ,(/ 'lev' /), 'I', ' ' ,'EXPLICITa' , gridname='geo_grid') ! part a
138 0 : call addfld ('EXPLICITb' ,(/ 'lev' /), 'I', ' ' ,'EXPLICITb' , gridname='geo_grid') ! part b
139 0 : call addfld ('EXPLICIT1' ,(/ 'lev' /), 'I', ' ' ,'EXPLICIT1' , gridname='geo_grid') ! complete
140 0 : call addfld ('EXPLICIT' ,(/ 'lev' /), 'I', ' ' ,'EXPLICIT' , gridname='geo_grid') ! final w/ poles
141 :
142 0 : call addfld ('EXPLICIT2' ,(/ 'lev' /), 'I', ' ' ,'EXPLICIT2' , gridname='geo_grid')
143 0 : call addfld ('EXPLICIT3' ,(/ 'lev' /), 'I', ' ' ,'EXPLICIT3' , gridname='geo_grid')
144 0 : call addfld ('TPHDZ0' ,(/ 'lev' /), 'I', ' ' ,'TPHDZ0' , gridname='geo_grid')
145 0 : call addfld ('TPHDZ1' ,(/ 'lev' /), 'I', ' ' ,'TPHDZ1' , gridname='geo_grid')
146 0 : call addfld ('DIVBZ' ,(/ 'lev' /), 'I', ' ' ,'DIVBZ' , gridname='geo_grid')
147 0 : call addfld ('HDZMBZ' ,(/ 'lev' /), 'I', ' ' ,'HDZMBZ' , gridname='geo_grid')
148 0 : call addfld ('HDZPBZ' ,(/ 'lev' /), 'I', ' ' ,'HDZPBZ' , gridname='geo_grid')
149 0 : call addfld ('P_COEFF0' ,(/ 'lev' /), 'I', ' ' ,'P_COEFF0' , gridname='geo_grid')
150 0 : call addfld ('Q_COEFF0' ,(/ 'lev' /), 'I', ' ' ,'Q_COEFF0' , gridname='geo_grid')
151 0 : call addfld ('R_COEFF0' ,(/ 'lev' /), 'I', ' ' ,'R_COEFF0' , gridname='geo_grid')
152 0 : call addfld ('P_COEFF0a' ,(/ 'lev' /), 'I', ' ' ,'P_COEFF0a' , gridname='geo_grid')
153 0 : call addfld ('Q_COEFF0a' ,(/ 'lev' /), 'I', ' ' ,'Q_COEFF0a' , gridname='geo_grid')
154 0 : call addfld ('DJINT' ,(/ 'lev' /), 'I', ' ' ,'DJINT' , gridname='geo_grid')
155 0 : call addfld ('BDOTU' ,(/ 'lev' /), 'I', ' ' ,'BDOTU' , gridname='geo_grid')
156 0 : call addfld ('R_COEFF0a' ,(/ 'lev' /), 'I', ' ' ,'R_COEFF0a' , gridname='geo_grid')
157 0 : call addfld ('P_COEFF1' ,(/ 'lev' /), 'I', ' ' ,'P_COEFF1' , gridname='geo_grid')
158 0 : call addfld ('Q_COEFF1' ,(/ 'lev' /), 'I', ' ' ,'Q_COEFF1' , gridname='geo_grid')
159 0 : call addfld ('R_COEFF1' ,(/ 'lev' /), 'I', ' ' ,'R_COEFF1' , gridname='geo_grid')
160 0 : call addfld ('P_COEFF2' ,(/ 'lev' /), 'I', ' ' ,'P_COEFF2' , gridname='geo_grid')
161 0 : call addfld ('Q_COEFF2' ,(/ 'lev' /), 'I', ' ' ,'Q_COEFF2' , gridname='geo_grid')
162 0 : call addfld ('R_COEFF2' ,(/ 'lev' /), 'I', ' ' ,'R_COEFF2' , gridname='geo_grid')
163 :
164 0 : call addfld ('P_COEFF' ,(/ 'lev' /), 'I', ' ' ,'P_COEFF' , gridname='geo_grid') ! final w/ poles
165 0 : call addfld ('Q_COEFF' ,(/ 'lev' /), 'I', ' ' ,'Q_COEFF' , gridname='geo_grid') ! final w/ poles
166 0 : call addfld ('R_COEFF' ,(/ 'lev' /), 'I', ' ' ,'R_COEFF' , gridname='geo_grid') ! final w/ poles
167 :
168 0 : call addfld ('OP_SOLVE' ,(/ 'lev' /), 'I', ' ' ,'OP_SOLVE' , gridname='geo_grid')
169 0 : call addfld ('OP_OUT' ,(/ 'lev' /), 'I', 'cm^3' ,'OPLUS (oplus_xport output)', gridname='geo_grid')
170 0 : call addfld ('OPNM_OUT' ,(/ 'lev' /), 'I', 'cm^3' ,'OPNM_OUT' , gridname='geo_grid')
171 0 : call addfld ('BMOD2' ,(/ 'lev' /), 'I', ' ' ,'BMOD2' , gridname='geo_grid')
172 :
173 0 : call addfld ('OPLUS_FLUX', horiz_only , 'I', ' ','OPLUS_FLUX', gridname='geo_grid')
174 0 : call addfld ('OPLUS_DIVB', horiz_only , 'I', ' ','OPLUS_DIVB', gridname='geo_grid')
175 0 : call addfld ('OPLUS_BX' , horiz_only , 'I', ' ','OPLUS_BX' , gridname='geo_grid')
176 0 : call addfld ('OPLUS_BY' , horiz_only , 'I', ' ','OPLUS_BY' , gridname='geo_grid')
177 0 : call addfld ('OPLUS_BZ' , horiz_only , 'I', ' ','OPLUS_BZ' , gridname='geo_grid')
178 0 : call addfld ('OPLUS_BMAG', horiz_only , 'I', ' ','OPLUS_BMAG', gridname='geo_grid')
179 : endif
180 :
181 2304 : allocate(zp(nlev)) ! log pressure (as in TIEGCM)
182 1536 : allocate(expz(nlev)) ! exp(-zp)
183 768 : zp = nan
184 768 : expz = nan
185 768 : end subroutine oplus_init
186 :
187 : !-----------------------------------------------------------------------
188 36480 : subroutine oplus_xport(tn, te, ti, un, vn, om, zg, o2, o1, n2, op_in, &
189 36480 : opnm_in, mbar, ui, vi, wi, pmid, op_out, opnm_out, &
190 : i0, i1, j0, j1, nspltop, ispltop)
191 : !
192 : ! All input fields from dpie_coupling are in "TIEGCM" format, i.e.,
193 : ! longitude (-180->180), vertical (bot2top), and units (CGS).
194 : !
195 768 : use edyn_mpi, only: mp_geo_halos,mp_pole_halos,setpoles
196 : use edyn_geogrid, only: glat, nlat, nlev
197 : use trsolv_mod, only: trsolv
198 : !
199 : ! Transport O+ ion.
200 : ! March-May, 2015 B.Foster: Adapted from TIEGCM (oplus.F) for WACCM-X.
201 : !
202 : ! Notes:
203 : ! - waccmx_opt='ionosphere' must be set in user_nl_cam for te,ti inputs to have values
204 : !
205 : ! Args:
206 : !
207 : integer,intent(in) :: &
208 : i0, & ! grid%ifirstxy
209 : i1, & ! grid%ilastxy
210 : j0, & ! grid%jfirstxy
211 : j1 ! grid%jlastxy
212 : integer,intent(in) :: nspltop,ispltop
213 : !
214 : ! Input fields without halo points (lon +/-180, vertical bot2top, CGS units):
215 : !
216 : real(r8),intent(in) :: tn (nlev,i0-2:i1+2,j0-2:j1+2) ! neutral temperature (deg K)
217 : real(r8),intent(in) :: te (nlev,i0-2:i1+2,j0-2:j1+2) ! electron temperature (deg K)
218 : real(r8),intent(in) :: ti (nlev,i0-2:i1+2,j0-2:j1+2) ! ion temperature (deg K)
219 : real(r8),intent(in) :: un (nlev,i0-2:i1+2,j0-2:j1+2) ! neutral zonal wind (cm/s)
220 : real(r8),intent(in) :: vn (nlev,i0-2:i1+2,j0-2:j1+2) ! neutral meridional wind (cm/s)
221 : real(r8),intent(in) :: om (nlev,i0-2:i1+2,j0-2:j1+2) ! omega (1/s)
222 : real(r8),intent(in) :: o2 (nlev,i0-2:i1+2,j0-2:j1+2) ! o2 (mmr)
223 : real(r8),intent(in) :: o1 (nlev,i0-2:i1+2,j0-2:j1+2) ! o (mmr)
224 : real(r8),intent(in) :: n2 (nlev,i0-2:i1+2,j0-2:j1+2) ! n2 (mmr)
225 : real(r8),intent(in) :: mbar (nlev,i0-2:i1+2,j0-2:j1+2) ! mean molecular weight
226 :
227 : real(r8),intent(in) :: op_in(nlev,i0:i1,j0:j1) ! O+ density (cm^3)
228 : real(r8),intent(in) :: opnm_in(nlev,i0:i1,j0:j1) ! O+ density (cm^3) at time-1
229 : real(r8),intent(in) :: zg (nlev,i0:i1,j0:j1) ! geopotential height (cm)
230 : !
231 : ! Ion drifts from edynamo (also in tiegcm-format):
232 : !
233 : real(r8),intent(in) :: ui(nlev,i0:i1,j0:j1) ! zonal ion drift
234 : real(r8),intent(in) :: vi(nlev,i0:i1,j0:j1) ! meridional ion drift
235 : real(r8),intent(in) :: wi(nlev,i0:i1,j0:j1) ! vertical ion drift
236 : real(r8),intent(in) :: pmid(nlev) ! pressure at midpoints (Pa)
237 : !
238 : ! Output:
239 : !
240 : real(r8),intent(out) :: &
241 : op_out (nlev,i0:i1,j0:j1), & ! O+ output
242 : opnm_out(nlev,i0:i1,j0:j1) ! O+ output at time n-1
243 : !
244 : ! Local:
245 : !
246 : integer :: i, j, k, lat, jm1, jp1, jm2, jp2, lat0, lat1
247 : real(r8), dimension(i0:i1,j0:j1) :: &
248 72960 : opflux, & ! upward number flux of O+ (returned by sub oplus_flux)
249 72960 : dvb ! divergence of B-field
250 : !
251 : ! Local inputs with added halo points in lat,lon:
252 : !
253 72960 : real(r8),dimension(nlev,i0-2:i1+2,j0-2:j1+2),target :: op, opnm
254 :
255 : real(r8),dimension(nlev,i0-2:i1+2,j0-2:j1+2),target :: &
256 72960 : tr ,& ! Reduced temperature (.5*(tn+ti))
257 72960 : tp ,& ! Plasma temperature N(O+)*(te+ti)
258 72960 : dj ,& ! diffusion coefficients
259 72960 : bvel ,& ! bvel @ j = (B.U)*N(O+)
260 72960 : diffj ,& ! (D/(H*DZ)*2.*TP+M*G/R)*N(O+)
261 72960 : bdotdh_op ,& ! (b(h)*del(h))*phi
262 72960 : bdotdh_opj ,& ! (b(h)*del(h))*phi
263 72960 : bdotdh_diff ,& ! (b(h)*del(h))*phi
264 72960 : opnm_smooth ! O+ at time-1, smoothed
265 :
266 : real(r8), dimension(nlev,i0:i1,j0:j1) :: & ! for saving to histories
267 72960 : diag0, diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, &
268 72960 : diag10, diag11, diag12, diag13, diag14, diag15, diag16, diag17, &
269 72960 : diag18, diag19, diag20, diag21, diag22, diag23, diag24, diag25, &
270 72960 : diag26, diag27
271 72960 : real(r8), dimension(nlev,i0:i1,j0-1:j1+1) :: hj ! scale height
272 : real(r8) :: gmr,dtime,dtx2,dtx2inv
273 : real(r8), dimension(nlev,i0:i1) :: &
274 72960 : bdzdvb_op, &
275 72960 : tp1, &
276 72960 : divbz
277 : real(r8),dimension(nlev,i0:i1,j0:j1) :: &
278 72960 : hdz, &
279 72960 : tphdz0, &
280 72960 : tphdz1, &
281 72960 : djint, &
282 72960 : hdzmbz, &
283 72960 : hdzpbz, &
284 72960 : bdotu
285 : ! for term analysis, lei, 07
286 : real(r8),dimension(nlev,i0:i1,j0:j1) :: &
287 72960 : op_dt, & ! dn/dt
288 72960 : amb_diff,& ! ambipole diffion
289 72960 : dwind, & ! neutral wind transport
290 72960 : dfield ! electric field transport
291 : !
292 : ! Arguments for tridiagonal solver trsolv (no halos):
293 : real(r8), dimension(nlev,i0:i1,j0:j1) :: &
294 109440 : explicit, explicit_a, explicit_b, p_coeff, q_coeff, r_coeff
295 :
296 72960 : real(r8), dimension(i0:i1) :: ubca, ubcb ! O+ upper boundary
297 : real(r8), parameter :: one=1._r8
298 : logical :: calltrsolv
299 : !
300 : ! Pointers for multiple-field calls (e.g., mp_geo_halos)
301 : integer :: nfields
302 36480 : real(r8), allocatable :: polesign(:)
303 36480 : type(array_ptr_type), allocatable :: ptrs(:)
304 :
305 72960 : real(r8) :: zpmid(nlev), opfloor
306 : real(r8), parameter :: opmin=3000.0_r8
307 : !
308 : ! Execute:
309 : !
310 36480 : dtime = get_step_size() ! step size in seconds
311 36480 : dtime = dtime / dble(nspltop)
312 36480 : dtx2 = 2._r8*dtime
313 36480 : dtx2inv = 1._r8/dtx2
314 :
315 36480 : if ((is_first_step() .or. is_first_restart_step()) .and. ispltop==1) then
316 768 : if (masterproc) then
317 : write(iulog,"(a,es12.4,a,es12.4,a,es12.4)") &
318 2 : 'oplus: shapiro=', shapiro_const, ', dtsmooth=', dtsmooth, &
319 4 : ', dtsmooth_div2=', dtsmooth_div2
320 2 : write(iulog,"('oplus: shr_const_g=',f8.3)") shr_const_g
321 : end if
322 : end if
323 :
324 : !
325 : ! zp,expz are declared in edyn_geogrid.F90, and allocated in sub
326 : ! set_geogrid (edyn_init.F90). pmid was passed in here (bot2top)
327 : ! from dpie_coupling.
328 : !
329 : ! kbot is the k-index at the bottom of O+ transport calculations,
330 : ! corresponding to pressure pbot.
331 : !
332 36480 : if ((is_first_step().or.is_first_restart_step()).and.ispltop==1) then
333 65280 : kloop: do k=1,nlev
334 65280 : if ( pmid(k) <= pbot) then
335 768 : kbot = k
336 768 : exit kloop
337 : end if
338 : enddo kloop
339 100608 : do k=1,nlev
340 99840 : zp(k) = -log(pmid(k)*10._r8/p0)
341 100608 : expz(k) = exp(-zp(k))
342 : enddo
343 : if (debug.and.masterproc) then
344 : write(iulog,"(a,i4,a,es12.4,a,es12.4)") &
345 : 'oplus: kbot=', kbot, ', pmid(kbot)=', pmid(kbot), &
346 : ', zp(kbot)=', zp(kbot)
347 : end if
348 : end if
349 :
350 36480 : if (kbot < 1) then
351 0 : call endrun('oplus_xport: kbot is not set')
352 : endif
353 :
354 36480 : dzp = zp(nlev) - zp(nlev-1) ! use top 2 levels (typically dzp=0.5)
355 :
356 : if (debug .and. masterproc) then
357 : write(iulog,"('oplus: nlev=',i3,' zp (bot2top) =',/,(6es12.3))") nlev, zp
358 : write(iulog,"('oplus: nlev=',i3,' expz (bot2top) =',/,(6es12.3))") nlev, expz
359 : write(iulog,"('oplus: nlev=',i3,' dzp =',/,(6es12.3))") nlev, dzp
360 : end if
361 : !
362 : ! Set subdomain blocks from input (composition is in mmr):
363 : !
364 : !$omp parallel do private(i, j, k)
365 4778880 : do k = 1, nlev
366 19006080 : do j = j0, j1
367 189696000 : do i = i0, i1
368 170726400 : op(k,i,j) = op_in(k,i,j)
369 184953600 : opnm(k,i,j) = opnm_in(k,i,j)
370 : end do
371 : end do
372 : end do
373 :
374 : !
375 : ! Define halo points on inputs:
376 : ! WACCM has global longitude values at the poles (j=1,j=nlev)
377 : ! (they are constant for most, except the winds.)
378 : !
379 : ! Set two halo points in lat,lon:
380 : ! real(r8),dimension(nlev,i0-2:i1+2,j0-2:j1+2),target :: tn,te,etc.
381 : !
382 36480 : nfields = 2
383 36480 : allocate(ptrs(nfields),polesign(nfields))
384 :
385 36480 : ptrs(1)%ptr => op ; ptrs(2)%ptr => opnm
386 109440 : polesign = 1._r8
387 : !
388 : ! mp_geo_halos first arg:
389 : ! type(array_ptr_type) :: fmsub(nf) ! (lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2)
390 : !
391 36480 : call mp_geo_halos(ptrs,1,nlev,i0,i1,j0,j1,nfields)
392 : !
393 : ! Set latitude halo points over the poles (this does not change the poles).
394 : ! (the 2nd halo over the poles will not actually be used (assuming lat loops
395 : ! are lat=2,nlat-1), because jp1,jm1 will be the pole itself, and jp2,jm2
396 : ! will be the first halo over the pole)
397 : !
398 : ! mp_pole_halos first arg:
399 : ! type(array_ptr_type) :: f(nf) ! (nlev,i0-2:i1+2,j0-2:j1+2)
400 :
401 36480 : call mp_pole_halos(ptrs,1,nlev,i0,i1,j0,j1,nfields,polesign)
402 36480 : deallocate(ptrs,polesign)
403 :
404 : !
405 : ! Use below to exclude the poles (lat=2,nlat-1) from latitude scans.
406 : !
407 36480 : lat0 = j0
408 36480 : lat1 = j1
409 36480 : if (j0 == 1) lat0 = 2
410 36480 : if (j1 == nlat) lat1 = nlat-1
411 : !
412 : ! Save input fields to WACCM histories. Sub savefld_waccm_switch converts
413 : ! fields from tiegcm-format to waccm-format before saving to waccm histories.
414 : !
415 36480 : if (debug_hist) then
416 0 : call savefld_waccm(tn(:,i0:i1,j0:j1),'OPLUS_TN',nlev,i0,i1,j0,j1)
417 0 : call savefld_waccm(te(:,i0:i1,j0:j1),'OPLUS_TE',nlev,i0,i1,j0,j1)
418 0 : call savefld_waccm(ti(:,i0:i1,j0:j1),'OPLUS_TI',nlev,i0,i1,j0,j1)
419 0 : call savefld_waccm(un(:,i0:i1,j0:j1),'OPLUS_UN',nlev,i0,i1,j0,j1)
420 0 : call savefld_waccm(vn(:,i0:i1,j0:j1),'OPLUS_VN',nlev,i0,i1,j0,j1)
421 0 : call savefld_waccm(om(:,i0:i1,j0:j1),'OPLUS_OM',nlev,i0,i1,j0,j1)
422 0 : call savefld_waccm(zg(:,i0:i1,j0:j1),'OPLUS_Z' ,nlev,i0,i1,j0,j1)
423 0 : call savefld_waccm(o2(:,i0:i1,j0:j1),'OPLUS_O2',nlev,i0,i1,j0,j1)
424 0 : call savefld_waccm(o1(:,i0:i1,j0:j1),'OPLUS_O1',nlev,i0,i1,j0,j1)
425 0 : call savefld_waccm(n2(:,i0:i1,j0:j1),'OPLUS_N2',nlev,i0,i1,j0,j1)
426 0 : call savefld_waccm(op(:,i0:i1,j0:j1),'OPLUS_OP',nlev,i0,i1,j0,j1)
427 0 : call savefld_waccm(ui(:,i0:i1,j0:j1),'OPLUS_UI',nlev,i0,i1,j0,j1)
428 0 : call savefld_waccm(vi(:,i0:i1,j0:j1),'OPLUS_VI',nlev,i0,i1,j0,j1)
429 0 : call savefld_waccm(wi(:,i0:i1,j0:j1),'OPLUS_WI',nlev,i0,i1,j0,j1)
430 0 : call savefld_waccm(mbar(:,i0:i1,j0:j1),'OPLUS_MBAR',nlev,i0,i1,j0,j1)
431 0 : call savefld_waccm(opnm(:,i0:i1,j0:j1),'OPLUS_OPNM',nlev,i0,i1,j0,j1)
432 : endif
433 : !
434 : ! Initialize output op_out with input op at 1:kbot-1, to retain values from
435 : ! bottom of column up to kbot. This routine will change (transport) these
436 : ! outputs only from kbot to the top (nlev).
437 : !
438 172185600 : op_out = 0._r8
439 172185600 : opnm_out = 0._r8
440 111774720 : op_out (1:kbot-1,i0:i1,j0:j1) = op (1:kbot-1,i0:i1,j0:j1)
441 111774720 : opnm_out(1:kbot-1,i0:i1,j0:j1) = opnm(1:kbot-1,i0:i1,j0:j1)
442 : !
443 : ! Sub oplus_flux returns upward number flux of O+ in opflux
444 : ! Output opflux(i,j) is 2d lon x lat subdomain:
445 : !
446 36480 : call oplus_flux(opflux,i0,i1,j0,j1)
447 36480 : if (debug_hist) then
448 0 : call savefld_waccm(opflux(i0:i1,j0:j1),'OPLUS_FLUX',1,i0,i1,j0,j1)
449 : endif
450 : !
451 : ! Divergence of B (mag field) is returned by divb in dvb(i0:i1,j0:j1)
452 : !
453 36480 : call divb(dvb,i0,i1,j0,j1)
454 36480 : if (debug_hist) then
455 0 : call savefld_waccm(dvb(i0:i1,j0:j1),'OPLUS_DIVB',1,i0,i1,j0,j1)
456 : endif
457 : !
458 : ! The solver will be called only if calltrsolv=true. It is sometimes
459 : ! set false when skipping parts of the code for debug purposes.
460 : !
461 : calltrsolv = .true.
462 :
463 535526400 : tr = 0._r8
464 535526400 : tp = 0._r8
465 535526400 : dj = 0._r8
466 286951680 : hj = 0._r8
467 535526400 : bvel= 0._r8
468 535526400 : diffj = 0._r8
469 535526400 : opnm_smooth = 0._r8
470 172185600 : diag0 =0._r8
471 36480 : grav_cm = shr_const_g * 100._r8 ! m/s^2 -> cm/s^2
472 : !
473 : !----------------------- Begin first latitude scan ---------------------
474 143640 : do lat=lat0,lat1
475 107160 : jm2 = lat-2
476 107160 : jm1 = lat-1
477 107160 : jp1 = lat+1
478 107160 : jp2 = lat+2
479 : !
480 : ! as of April, 2015, TIEGCM incorrectly uses te+ti instead of tn+ti
481 : ! This has not been fixed in TIEGCM, because fixing it causes a tuning
482 : ! problem (ask Hanli and Wenbin). For WACCM, it is correct as below.
483 : ! (see also tp)
484 : !
485 : !$omp parallel do private(i,k)
486 1393080 : do i=i0,i1
487 : !
488 : ! Reduced temperature (tpj in tiegcm):
489 : ! 'OPLUS_TR' (has constants at poles)
490 : !
491 60545400 : do k=kbot,nlev
492 59152320 : tr(k,i,jm1) = 0.5_r8*(tn(k,i,jm1)+ti(k,i,jm1))
493 59152320 : tr(k,i,lat) = 0.5_r8*(tn(k,i,lat)+ti(k,i,lat))
494 60438240 : tr(k,i,jp1) = 0.5_r8*(tn(k,i,jp1)+ti(k,i,jp1))
495 : enddo
496 : enddo ! i=i0,i1
497 : !
498 : ! rrk returns ambipolar diffusion coefficients in d(jm1),dj(lat),djp1(jp1):
499 : ! 'OPLUS_DJ' (has constants at poles)
500 : !
501 : call rrk( &
502 214320 : tn(kbot:nlev,i0:i1,jm1),mbar(kbot:nlev,i0:i1,jm1), &
503 214320 : o2(kbot:nlev,i0:i1,jm1),o1 (kbot:nlev,i0:i1,jm1), &
504 214320 : n2(kbot:nlev,i0:i1,jm1),tr (kbot:nlev,i0:i1,jm1), &
505 423817800 : dj(kbot:nlev,i0:i1,jm1),i0,i1,kbot,nlev)
506 :
507 : call rrk( &
508 214320 : tn(kbot:nlev,i0:i1,lat),mbar(kbot:nlev,i0:i1,lat), &
509 214320 : o2(kbot:nlev,i0:i1,lat),o1 (kbot:nlev,i0:i1,lat), &
510 214320 : n2(kbot:nlev,i0:i1,lat),tr (kbot:nlev,i0:i1,lat), &
511 423817800 : dj(kbot:nlev,i0:i1,lat),i0,i1,kbot,nlev)
512 :
513 : call rrk( &
514 214320 : tn(kbot:nlev,i0:i1,jp1),mbar(kbot:nlev,i0:i1,jp1), &
515 214320 : o2(kbot:nlev,i0:i1,jp1),o1 (kbot:nlev,i0:i1,jp1), &
516 214320 : n2(kbot:nlev,i0:i1,jp1),tr (kbot:nlev,i0:i1,jp1), &
517 423817800 : dj(kbot:nlev,i0:i1,jp1),i0,i1,kbot,nlev)
518 : !
519 : ! Plasma temperature:
520 : ! 'OPLUS_TP0' (tp will get poles from jm1 and jp1)
521 : !
522 : !$omp parallel do private(i,k)
523 1393080 : do i=i0,i1
524 60545400 : do k=kbot,nlev
525 59152320 : tp(k,i,jm1) = te(k,i,jm1)+ti(k,i,jm1)
526 59152320 : tp(k,i,lat) = te(k,i,lat)+ti(k,i,lat)
527 60438240 : tp(k,i,jp1) = te(k,i,jp1)+ti(k,i,jp1)
528 : enddo
529 : enddo
530 60545400 : diag0(kbot:nlev,i0:i1,lat) = tp(kbot:nlev,i0:i1,lat)
531 : !
532 : ! Add poles to diag0:
533 750120 : if (j0==1.and.lat==2) diag0(kbot:nlev,i0:i1,j0) = tp(kbot:nlev,i0:i1,jm1)
534 750120 : if (j1==nlat.and.lat==nlat-1) diag0(kbot:nlev,i0:i1,j1) = tp(kbot:nlev,i0:i1,jp1)
535 : !
536 : ! Neutral scale height:
537 : ! 'OPLUS_HJ' (has constants at poles)
538 : !
539 : !$omp parallel do private(i,k)
540 1393080 : do i=i0,i1
541 60545400 : do k=kbot,nlev
542 59152320 : hj(k,i,jm1) = gask * tn(k,i,jm1) / (mbar(k,i,jm1) * grav_cm)
543 59152320 : hj(k,i,lat) = gask * tn(k,i,lat) / (mbar(k,i,lat) * grav_cm)
544 60438240 : hj(k,i,jp1) = gask * tn(k,i,jp1) / (mbar(k,i,jp1) * grav_cm)
545 : enddo
546 : enddo
547 : !
548 : ! bvel @ jm1 = (B.U)*N(O+) (J-1)
549 : ! bvel @ j = (B.U)*N(O+) (J)
550 : ! bvel @ jp1 = (B.U)*N(O+) (J+1)
551 : ! 'OPLUS_BVEL' (has constants at poles)
552 : !
553 : ! Note bx,by,bz were set globally for all tasks by sub magfield
554 : ! (getapex.F90)
555 : !
556 : !$omp parallel do private(i,k)
557 1393080 : do i=i0,i1
558 60545400 : do k=kbot,nlev
559 177456960 : bvel(k,i,jm1) = &
560 0 : (bx(i,jm1)*un(k,i,jm1)+by(i,jm1)*vn(k,i,jm1)+ &
561 236609280 : hj(k,i,jm1)*bz(i,jm1)*om(k,i,jm1))*op(k,i,jm1)
562 : bvel(k,i,lat) = &
563 0 : (bx(i,lat)*un(k,i,lat)+by(i,lat)*vn(k,i,lat)+ &
564 59152320 : hj(k,i,lat)*bz(i,lat)*om(k,i,lat))*op(k,i,lat)
565 59152320 : bvel(k,i,jp1) = &
566 0 : (bx(i,jp1)*un(k,i,jp1)+by(i,jp1)*vn(k,i,jp1)+ &
567 119590560 : hj(k,i,jp1)*bz(i,jp1)*om(k,i,jp1))*op(k,i,jp1)
568 : enddo ! k=kbot,nlev
569 : enddo ! i=lon0,lon1
570 : !
571 : ! Ambipolar diffusion is returned in diffj:
572 : ! 'OPLUS_DIFFJ' (will have constants at poles after this lat scan)
573 : !
574 321480 : call diffus(tp(kbot:nlev,i0:i1,jm1),op(kbot:nlev,i0:i1,jm1),hj(kbot:nlev,:,jm1), &
575 242181600 : diffj(kbot:nlev,i0:i1,jm1),i0,i1,kbot,nlev,lat)
576 321480 : call diffus(tp(kbot:nlev,i0:i1,lat),op(kbot:nlev,i0:i1,lat),hj(kbot:nlev,:,lat), &
577 242181600 : diffj(kbot:nlev,i0:i1,lat),i0,i1,kbot,nlev,lat)
578 321480 : call diffus(tp(kbot:nlev,i0:i1,jp1),op(kbot:nlev,i0:i1,jp1),hj(kbot:nlev,:,jp1), &
579 242181600 : diffj(kbot:nlev,i0:i1,jp1),i0,i1,kbot,nlev,lat)
580 : !
581 : ! 'OPLUS_TP1' (constants at the poles)
582 : !
583 : !$omp parallel do private(i,k)
584 1393080 : do i=i0,i1
585 60545400 : do k=kbot,nlev
586 59152320 : tp(k,i,jm2) = op(k,i,jm2)*(te(k,i,jm2)+ti(k,i,jm2))
587 59152320 : tp(k,i,jm1) = tp(k,i,jm1)*op(k,i,jm1)
588 59152320 : tp(k,i,lat) = tp(k,i,lat)*op(k,i,lat)
589 59152320 : tp(k,i,jp1) = tp(k,i,jp1)*op(k,i,jp1)
590 60438240 : tp(k,i,jp2) = op(k,i,jp2)*(te(k,i,jp2)+ti(k,i,jp2))
591 : enddo
592 : enddo
593 : !
594 : ! Latidinal shapiro smoother: opnm is O+ at time n-1.
595 : ! opnm_smooth will be used in explicit terms below.
596 : ! Smooth in latitude:
597 : ! 'OPNM_SMOOTH' (zero at poles)
598 : !
599 : !$omp parallel do private(i,k)
600 1429560 : do i=i0,i1
601 60545400 : do k=kbot,nlev
602 177456960 : opnm_smooth(k,i,lat) = opnm(k,i,lat)-shapiro_const* &
603 118304640 : (opnm(k,i,jp2)+opnm(k,i,jm2)-4._r8* &
604 118304640 : (opnm(k,i,jp1)+opnm(k,i,jm1))+6._r8* &
605 474504480 : opnm(k,i,lat))
606 : enddo ! k=kbot,nlev
607 : enddo ! i=i0,i1
608 : enddo ! end first latitude scan (lat=lat0,lat1)
609 : !
610 : !------------------------- End first latitude scan ---------------------
611 : !
612 : ! Set pole values for opnm_smooth. Do this before savefld calls, so plots will
613 : ! include the poles. All other fields in 1st lat scan got values at the poles
614 : ! via jm1,jp1 above.
615 : !
616 123703680 : call setpoles(opnm_smooth(kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
617 : !
618 : ! Save to history file (exclude halo points)
619 : !
620 36480 : if (debug_hist) then
621 0 : call savefld_waccm(tr (:,i0:i1,j0:j1),'OPLUS_TR' ,nlev,i0,i1,j0,j1)
622 0 : call savefld_waccm(dj (:,i0:i1,j0:j1),'OPLUS_DJ' ,nlev,i0,i1,j0,j1)
623 0 : call savefld_waccm(hj (:,i0:i1,j0:j1),'OPLUS_HJ' ,nlev,i0,i1,j0,j1)
624 0 : call savefld_waccm(bvel (:,i0:i1,j0:j1),'OPLUS_BVEL' ,nlev,i0,i1,j0,j1)
625 0 : call savefld_waccm(diffj(:,i0:i1,j0:j1),'OPLUS_DIFFJ',nlev,i0,i1,j0,j1)
626 0 : call savefld_waccm(diag0(:,i0:i1,j0:j1),'OPLUS_TP0' ,nlev,i0,i1,j0,j1)
627 0 : call savefld_waccm(tp (:,i0:i1,j0:j1),'OPLUS_TP1' ,nlev,i0,i1,j0,j1)
628 0 : call savefld_waccm(opnm_smooth(:,i0:i1,j0:j1),'OPNM_SMOOTH',nlev,i0,i1,j0,j1)
629 : endif
630 : !
631 : ! Set halo points where needed.
632 : !
633 36480 : nfields = 5
634 36480 : allocate(ptrs(nfields),polesign(nfields))
635 36480 : ptrs(1)%ptr => dj ; ptrs(2)%ptr => bvel ; ptrs(3)%ptr => diffj
636 36480 : ptrs(4)%ptr => tp ; ptrs(5)%ptr => opnm_smooth
637 218880 : polesign = 1._r8
638 :
639 36480 : call mp_geo_halos (ptrs,1,nlev,i0,i1,j0,j1,5)
640 36480 : call mp_pole_halos(ptrs,1,nlev,i0,i1,j0,j1,5,polesign)
641 :
642 36480 : deallocate(ptrs,polesign)
643 :
644 : !----------------------- Begin second latitude scan --------------------
645 535526400 : bdotdh_op = 0._r8
646 535526400 : bdotdh_opj = 0._r8
647 :
648 143640 : do lat=lat0,lat1
649 107160 : jm2 = lat-2
650 107160 : jm1 = lat-1
651 107160 : jp1 = lat+1
652 107160 : jp2 = lat+2
653 : !
654 : ! bdotdh_op = (B(H).DEL(H))*(D/(H*DZ)*TP+M*G/R)*N(O+)
655 : ! then bdotdh_op = d*bz*bdotdh_op
656 : ! real(r8),dimension(nlev,i0-2:i1+2,j0-2:j1+2) :: diffj
657 : ! real(r8),dimension(nlev,i0-2:i1+2,j0-2:j1+2) :: bdotdh_op
658 : ! 'BDOTDH_OP' (zero at the poles)
659 : !
660 : call bdotdh( &
661 107160 : diffj(kbot:nlev,i0:i1,jm1), &
662 214320 : diffj(kbot:nlev,:,lat ), & ! includes longitude halos
663 107160 : diffj(kbot:nlev,i0:i1,jp1), &
664 343019160 : bdotdh_op(kbot:nlev,i0:i1,lat),i0,i1,kbot,nlev,lat)
665 : !
666 : !$omp parallel do private( i, k )
667 1393080 : do i=i0,i1
668 60545400 : do k=kbot,nlev
669 60438240 : bdotdh_op(k,i,lat) = dj(k,i,lat)*bz(i,lat)*bdotdh_op(k,i,lat) ! BDOTDH_OP
670 : enddo ! k=kbot,nlev
671 : enddo ! i=i0,i1
672 : !
673 : ! bdotdh_opjm1 = (B(H).DEL(H))*2.*TP*N(O+) (J-1)
674 : ! bdotdh_opj = (B(H).DEL(H))*2.*TP*N(O+) (J)
675 : ! bdotdh_opjp1 = (B(H).DEL(H))*2.*TP*N(O+) (J+1)
676 : ! 'BDOTDH_OPJ' (has reasonable non-constant values at poles)
677 : !
678 : call bdotdh( &
679 107160 : tp(kbot:nlev,i0:i1,jm2), &
680 214320 : tp(kbot:nlev,:,jm1), &
681 107160 : tp(kbot:nlev,i0:i1,lat), &
682 343019160 : bdotdh_opj(kbot:nlev,i0:i1,jm1),i0,i1,kbot,nlev,jm1)
683 : call bdotdh( &
684 107160 : tp(kbot:nlev,i0:i1,jm1), &
685 214320 : tp(kbot:nlev,:,lat), &
686 107160 : tp(kbot:nlev,i0:i1,jp1), &
687 343019160 : bdotdh_opj(kbot:nlev,i0:i1,lat),i0,i1,kbot,nlev,lat)
688 : call bdotdh( &
689 107160 : tp(kbot:nlev,i0:i1,lat), &
690 214320 : tp(kbot:nlev,:,jp1), &
691 107160 : tp(kbot:nlev,i0:i1,jp2), &
692 343019160 : bdotdh_opj(kbot:nlev,i0:i1,jp1),i0,i1,kbot,nlev,jp1)
693 : !
694 : !$omp parallel do private( i, k )
695 1429560 : do i=i0,i1
696 60545400 : do k=kbot,nlev
697 59152320 : bdotdh_opj(k,i,jm1) = bdotdh_opj(k,i,jm1)*dj(k,i,jm1)
698 59152320 : bdotdh_opj(k,i,lat) = bdotdh_opj(k,i,lat)*dj(k,i,lat)
699 60438240 : bdotdh_opj(k,i,jp1) = bdotdh_opj(k,i,jp1)*dj(k,i,jp1)
700 : enddo ! k=kbot,nlev
701 : enddo ! i=i0,i1
702 : enddo ! lat=j0,j1 (end second lat scan)
703 : !
704 : !------------------------ End second latitude scan ---------------------
705 : !
706 : ! bdotdh_opj already has non-constant polar values, but bdotdh_op poles are zero.
707 : ! Sub setpoles will set poles to the zonal average of the latitude below each pole.
708 : !
709 : ! This may not be necessary, but do it for plotting:
710 123703680 : call setpoles(bdotdh_op(kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
711 :
712 36480 : if (debug_hist) then
713 0 : call savefld_waccm(bdotdh_op (:,i0:i1,j0:j1),'BDOTDH_OP' ,nlev,i0,i1,j0,j1)
714 0 : call savefld_waccm(bdotdh_opj(:,i0:i1,j0:j1),'BDOTDH_OPJ',nlev,i0,i1,j0,j1)
715 : endif
716 : !
717 : ! Note mp_geo_halos will overwrite jm1,jp1 that was set above.
718 : ! bdotdh_opj needs longitude halos for the bdotdh call below.
719 : !
720 : ! real(r8),dimension(nlev,i0-2:i1+2,j0-2:j1+2),target :: bdotdh_op,opj
721 : !
722 36480 : allocate(ptrs(1))
723 36480 : ptrs(1)%ptr => bdotdh_opj
724 36480 : call mp_geo_halos (ptrs,1,nlev,i0,i1,j0,j1,1)
725 36480 : call mp_pole_halos(ptrs,1,nlev,i0,i1,j0,j1,1,(/1._r8/))
726 36480 : deallocate(ptrs)
727 : !
728 : !----------------------- Begin third latitude scan ---------------------
729 : !
730 535526400 : bdotdh_diff = 0._r8
731 57383040 : bdzdvb_op = 0._r8
732 516483840 : explicit(1:nlev,i0:i1,j0:j1) = 0._r8 ; explicit_a(1:nlev,i0:i1,j0:j1)=0._r8 ; explicit_b(1:nlev,i0:i1,j0:j1)=0._r8
733 172185600 : hdz = 0._r8
734 172185600 : tphdz0 = 0._r8
735 172185600 : tphdz1 = 0._r8
736 172185600 : djint = 0._r8
737 57383040 : divbz = 0._r8
738 172185600 : hdzmbz = 0._r8
739 172185600 : hdzpbz = 0._r8
740 172185600 : p_coeff(1:nlev,i0:i1,j0:j1) = 0._r8
741 172185600 : q_coeff(1:nlev,i0:i1,j0:j1) = 0._r8
742 172185600 : r_coeff(1:nlev,i0:i1,j0:j1) = 0._r8
743 172185600 : bdotu = 0._r8
744 172185600 : op_dt = 0._r8
745 172185600 : amb_diff = 0._r8
746 172185600 : dwind = 0._r8
747 172185600 : dfield = 0._r8
748 :
749 860782080 : diag1 = 0._r8 ; diag2 = 0._r8 ; diag3 = 0._r8 ; diag4 = 0._r8 ; diag5 = 0._r8
750 860782080 : diag6 = 0._r8 ; diag7 = 0._r8 ; diag8 = 0._r8 ; diag9 = 0._r8 ; diag10= 0._r8
751 860782080 : diag11 = 0._r8 ; diag12= 0._r8 ; diag13= 0._r8 ; diag14= 0._r8 ; diag15= 0._r8
752 860782080 : diag16 = 0._r8 ; diag17= 0._r8 ; diag18= 0._r8 ; diag19= 0._r8 ; diag20= 0._r8
753 860782080 : diag21 = 0._r8 ; diag22= 0._r8 ; diag23= 0._r8 ; diag24= 0._r8 ; diag25= 0._r8
754 344334720 : diag26 = 0._r8 ; diag27= 0._r8
755 :
756 : !
757 : ! gmr = G*M(O+)/(2.*R)
758 : !
759 36480 : gmr = grav_cm*rmass_op/(2._r8*gask)
760 :
761 : !
762 : ! Globally, this loop is lat=2,nlat-1 (i.e., skipping the poles)
763 : !
764 143640 : do lat=lat0,lat1
765 107160 : jm2 = lat-2
766 107160 : jm1 = lat-1 ! this will be south pole for southern pes (j==1)
767 107160 : jp1 = lat+1 ! this will be north pole for northern pes (j==nlat)
768 107160 : jp2 = lat+2
769 : !
770 : ! bdotdh_opj = (B(H).DEL(H))*D*(B(H).DEL(H))*2.*TP*N(O+) (J)
771 : ! 'BDOTDH_DIFF' (zero at the poles)
772 : !
773 : call bdotdh( &
774 107160 : bdotdh_opj(kbot:nlev,i0:i1,jm1), &
775 214320 : bdotdh_opj(kbot:nlev,:,lat), & ! includes longitude halos
776 107160 : bdotdh_opj(kbot:nlev,i0:i1,jp1), &
777 343019160 : bdotdh_diff(kbot:nlev,i0:i1,lat),i0,i1,kbot,nlev,lat) ! BDOTDH_DIFF
778 : !
779 : ! bdzdvb_op = (BZ*D/(H*DZ)+DIV(*B))*S2
780 : ! bdzdvb returns bdzdvb_op(k,i).
781 : ! 'BDZDVB_OP' (zero at the poles)
782 : !
783 : ! real(r8),dimension(i0:i1,j0:j1) :: dvb
784 : ! real(r8),dimension(nlev,i0:i1,j0-1:j1+1) :: hj ! scale height
785 : ! real(r8),dimension(nlev,i0:i1) :: bdzdvb_op
786 :
787 : ! subroutine bdzdvb(phi,dvb,h,ans,lev0,lev1,lon0,lon1,lat)
788 : ! real(r8),intent(in) :: dvb(lon0:lon1)
789 : ! real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: phi,h
790 : ! real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: ans
791 : !
792 214320 : call bdzdvb(bdotdh_opj(kbot:nlev,i0:i1,lat),dvb(:,lat),hj(kbot:nlev,i0:i1,lat), &
793 181636200 : bdzdvb_op(kbot:nlev,i0:i1),kbot,nlev,i0,i1,lat)
794 168562680 : diag1(:,i0:i1,lat) = bdzdvb_op(:,i0:i1) ! BDZDVB_OP
795 : !
796 : ! Collect explicit terms:
797 : ! 'EXPLICIT0' (this will have poles set after third lat scan, before
798 : ! plotting. The poles will be constant in longitude, and
799 : ! may differ structurally from adjacent latitudes.
800 : !
801 : !$omp parallel do private( i, k )
802 1393080 : do i=i0,i1
803 60545400 : do k=kbot,nlev
804 59152320 : explicit(k,i,lat) = -one*(bdzdvb_op(k,i)+bdotdh_diff(k,i,lat)+bdotdh_op(k,i,lat))
805 60438240 : amb_diff(k,i,lat) = -explicit(k,i,lat)
806 : enddo ! k=kbot,nlev
807 : enddo ! i=i0,i1
808 168562680 : diag2(:,i0:i1,lat) = explicit(:,i0:i1,lat) ! EXPLICIT0
809 : !
810 : ! Ion drifts are interpolated to midpoints (is this necessary in WACCM?).
811 : !
812 : ! Need lon,lat halos for op, bvel, and bmod2
813 : ! op,bvel halos were set above, bmod2 was set in magfield (getapex.F90)
814 : ! (ui,vi,wi halos are not used here.)
815 : !
816 : ! bmod2 halos are set in sub magfield (getapex.F90), including nlat-1,nlat,nlat+1,
817 : ! and 1 halo point in longitude. Note bmod2 is global in lon and lat for all pe's.
818 : ! use getapex,only: bmod2 ! (0:nlonp1,jspole-1:jnpole+1)
819 : !
820 : ! When looping lat=2,nlat-1, this explicit has zero pole values,
821 : ! but there are still problems at processor longitude boundaries,
822 : ! especially near the south pole:
823 : ! 'EXPLICIT1' (zero at the poles)
824 : !
825 : !$omp parallel do private( i, k )
826 1393080 : do i=i0,i1
827 59259480 : do k=kbot,nlev-1
828 : !
829 : ! Original TIEGCM statement:
830 : ! explicit(k,i) = explicit(k,i)+1._r8/(2._r8*Rearth)* &
831 : ! (1._r8/(cs(lat)*dlamda)*(bx(i,lat)* &
832 : ! (bvel(k,i+1,lat)-bvel(k,i-1,lat))+ &
833 : ! 0.5_r8*(ui(k,i,lat)+ui(k+1,i,lat))*bmod2(i,lat)**2* &
834 : ! (op(k,i+1,lat)/bmod2(i+1,lat)**2- &
835 : ! op(k,i-1,lat)/bmod2(i-1,lat)**2))+ &
836 : !
837 : ! 1._r8/dphi*(by(i,lat)*(bvel(k,i,jp1)-bvel(k,i,jm1))+ &
838 : ! 0.5_r8*(vi(k,i,lat)+vi(k+1,i,lat))*bmod2(i,lat)**2* &
839 : ! (op(k,i,jp1)/bmod2(i,jp1)**2- &
840 : ! op(k,i,jm1)/bmod2(i,jm1)**2)))
841 : !
842 : ! Break it into two pieces and put together for debug:
843 : !
844 : ! 'EXPLICITa'
845 115732800 : explicit_a(k,i,lat) = (bx(i,lat)* &
846 173599200 : (bvel(k,i+1,lat)-bvel(k,i-1,lat))+ &
847 57866400 : 0.5_r8*(ui(k,i,lat)+ui(k+1,i,lat))*bmod2(i,lat)**2* &
848 57866400 : (op(k,i+1,lat)/bmod2(i+1,lat)**2- &
849 462931200 : op(k,i-1,lat)/bmod2(i-1,lat)**2))
850 : !
851 : ! 'EXPLICITb'
852 : !
853 : explicit_b(k,i,lat) = &
854 115732800 : (by(i,lat)*(bvel(k,i,jp1)-bvel(k,i,jm1))+ &
855 0 : 0.5_r8*(vi(k,i,lat)+vi(k+1,i,lat))*bmod2(i,lat)**2* &
856 57866400 : (op(k,i,jp1)/bmod2(i,jp1)**2- &
857 231465600 : op(k,i,jm1)/bmod2(i,jm1)**2))
858 : !
859 : ! 'EXPLICIT1'
860 : ! explicit will receive polar values after this latitude scan.
861 : !
862 : explicit(k,i,lat) = explicit(k,i,lat)+1._r8/(2._r8*Rearth)* &
863 0 : (1._r8/(cs(lat)*dlamda)*explicit_a(k,i,lat)+ &
864 57866400 : 1._r8/dphi*explicit_b(k,i,lat))
865 :
866 57866400 : dfield(k,i,lat) = -(explicit(k,i,lat)+amb_diff(k,i,lat))
867 : !
868 : ! explicit is bad at i=1,72,73,144 near south pole (npole appears to be ok)
869 : ! This does not appear to adversely affect the final O+ output, and TIEGCM
870 : ! has the same high magnitudes, so am ignoring this for now. The high magnitudes
871 : ! are near the south pole, at processor longitude boundaries (implicating an error
872 : ! with longitude halo points).
873 : !
874 1285920 : if (debug) then
875 : if (explicit(k,i,lat) < -300._r8 .or. explicit(k,i,lat) > 300._r8) then
876 : write(iulog,"('>>> bad explicit: k,i,lat=',3i4,' explicit=',es12.4)") &
877 : k,i,lat,explicit(k,i,lat)
878 : write(iulog,"(' cs(lat) =',3es12.4)") cs(lat)
879 : write(iulog,"(' op(k,i-1:i+1,lat) =',3es12.4)") op(k,i-1:i+1,lat)
880 : write(iulog,"(' op(k,i,jm1:jp1) =',3es12.4)") op(k,i,jm1:jp1)
881 : write(iulog,"(' bvel(k,i-1:i+1,lat)=',3es12.4)") bvel(k,i-1:i+1,lat)
882 : write(iulog,"(' bvel(k,i,jm1:jp1) =',3es12.4)") bvel(k,i,jm1:jp1)
883 : write(iulog,"(' bmod2(i-1:i+1,lat) =',3es12.4)") bmod2(i-1:i+1,lat)
884 : write(iulog,"(' bmod2(i,jm1:jp1) =',3es12.4)") bmod2(i,jm1:jp1)
885 : write(iulog,"(' ui(k:k+1,i,lat) =',2es12.4)") ui(k:k+1,i,lat)
886 : write(iulog,"(' vi(k:k+1,i,lat) =',2es12.4)") vi(k:k+1,i,lat)
887 : write(iulog,"(' bx,by(i,lat) =',2es12.4)") bx(i,lat),by(i,lat)
888 : endif
889 : endif
890 :
891 : enddo ! k=kbot,nlev-1
892 : enddo ! i=i0,i1
893 :
894 : !$omp parallel do private( k )
895 5036520 : do k=kbot,nlev
896 64188840 : diag25(k,i0:i1,lat) = bmod2(i0:i1,lat) ! BMOD2 (redundant in vertical)
897 : enddo
898 168562680 : diag26(:,i0:i1,lat) = explicit_a(:,i0:i1,lat) ! EXPLICITa
899 168562680 : diag27(:,i0:i1,lat) = explicit_b(:,i0:i1,lat) ! EXPLICITb
900 168562680 : diag3 (:,i0:i1,lat) = explicit (:,i0:i1,lat) ! EXPLICIT1
901 :
902 : !$omp parallel do private( i )
903 1393080 : do i=i0,i1
904 1393080 : dvb(i,lat) = dvb(i,lat)/bz(i,lat)
905 : enddo ! i=i0,i1
906 :
907 : !$omp parallel do private( i, k )
908 1393080 : do i=i0,i1
909 60545400 : do k=kbot,nlev
910 59152320 : hdz(k,i,lat) = 1._r8/(hj(k,i,lat)*dzp)
911 60438240 : tp1(k,i) = 0.5_r8*(ti(k,i,lat)+te(k,i,lat))
912 : enddo ! k=kbot,nlev
913 : enddo ! i=i0,i1
914 :
915 : !$omp parallel do private( i, k )
916 1393080 : do i=i0,i1
917 59259480 : do k=kbot,nlev-1
918 57866400 : tphdz1(k+1,i,lat) = 2._r8*tp1(k+1,i)*(0.5_r8*(hdz(k,i,lat)+hdz(k+1,i,lat)))+gmr
919 59152320 : tphdz0(k+1,i,lat) = 2._r8*tp1(k ,i)*(0.5_r8*(hdz(k,i,lat)+hdz(k+1,i,lat)))-gmr
920 : enddo ! k=kbot,nlev-1
921 : enddo ! i=lon0,lon1
922 : !
923 : ! Upper and lower boundaries:
924 : ! Both TPHDZ0 and TPHDZ1 are zero at the poles.
925 : !
926 : ! 5/9/15: Appears to be a problem in TPHDZ0,1 near kbot, maybe is zero?
927 : !
928 : !$omp parallel do private( i )
929 1393080 : do i=i0,i1
930 2571840 : tphdz1(kbot,i,lat) = 2._r8*tp1(kbot,i)* &
931 3857760 : (1.5_r8*hdz(kbot,i,lat)-0.5_r8*hdz(kbot+1,i,lat))+gmr
932 3857760 : tphdz1(nlev,i,lat) = 2._r8*(2._r8*tp1(nlev-1,i)-tp1(nlev-2,i))* &
933 3857760 : (1.5_r8*hdz(nlev-1,i,lat)-0.5_r8*hdz(nlev-2,i,lat))+gmr
934 : tphdz0(kbot,i,lat) = 2._r8*(2._r8*tp1(kbot,i)-tp1(kbot+1,i))* &
935 1285920 : (1.5_r8*hdz(kbot,i,lat)-0.5_r8*hdz(kbot+1,i,lat))-gmr
936 : tphdz0(nlev,i,lat) = 2._r8*tp1(nlev-1,i)* &
937 1393080 : (1.5_r8*hdz(nlev-1,i,lat)-0.5_r8*hdz(nlev-2,i,lat))-gmr
938 : enddo ! i=i0,i1
939 168562680 : diag4(:,i0:i1,lat) = tphdz0(:,i0:i1,lat) ! TPHDZ0
940 168562680 : diag5(:,i0:i1,lat) = tphdz1(:,i0:i1,lat) ! TPHDZ1
941 : !
942 : ! djint = dj diffusion at interfaces:
943 : ! 'DJINT' (zero at the poles - messes up the plots - may give
944 : ! diag6 polar values after the lat scan, before plotting)
945 : !
946 : !$omp parallel do private( i, k )
947 1393080 : do i=i0,i1
948 59152320 : do k=kbot,nlev-1
949 59152320 : djint(k+1,i,lat) = 0.5_r8*(dj(k,i,lat)+dj(k+1,i,lat))
950 : enddo
951 1285920 : djint(kbot,i,lat) = (1.5_r8*dj(kbot ,i,lat)-0.5_r8*dj(kbot+1,i,lat))
952 1393080 : djint(nlev,i,lat) = (1.5_r8*dj(nlev-1,i,lat)-0.5_r8*dj(nlev-2,i,lat))
953 : enddo ! i=i0,i1
954 168562680 : diag6(:,i0:i1,lat) = djint(:,i0:i1,lat) ! DJINT
955 : !
956 : ! divbz = (DIV(B)+(DH*D*BZ)/(D*BZ)
957 : ! 'DIVBZ' Field appears as a line following mins along magnetic equator (zero at poles)
958 : ! Field may be zero at kbot? Lat slices look strange between +/- 14 deg lat.
959 : !
960 : !$omp parallel do private( i, k )
961 1393080 : do i=i0,i1
962 60545400 : do k=kbot,nlev
963 118304640 : divbz(k,i) = &
964 59152320 : dvb(i,lat)+1._r8/(Rearth*dj(k,i,lat)*bz(i,lat)**2)*(bx(i,lat)/ &
965 118304640 : cs(lat)*(dj(k,i+1,lat)*bz(i+1,lat)-dj(k,i-1,lat)* &
966 118304640 : bz(i-1,lat))/(2._r8*dlamda)+by(i,lat)*(dj(k,i,jp1)* &
967 356199840 : bz(i,jp1)-dj(k,i,jm1)*bz(i,jm1))/(2._r8*dphi))
968 : enddo ! k=kbot,nlev
969 : enddo ! i=i0,i1
970 168562680 : diag7(:,i0:i1,lat) = divbz(:,i0:i1) ! DIVBZ
971 : !
972 : ! hdzmbz = (1./(H*DZ)-(DIV(B)+DH*D*BZ/(D*BZ))/(2*BZ))*BZ**2
973 : ! hdzpbz = (1./(H*DZ)+(DIV(B)+DH*D*BZ/(D*BZ))/(2*BZ))*BZ**2
974 : ! 'HDZMBZ' and 'HDZPBZ' are zero at the poles.
975 : !
976 : !$omp parallel do private( i, k )
977 1393080 : do i=i0,i1
978 60545400 : do k=kbot,nlev
979 59152320 : hdzmbz(k,i,lat) = (hdz(k,i,lat)-0.5_r8*divbz(k,i))*bz(i,lat)**2
980 60438240 : hdzpbz(k,i,lat) = (hdz(k,i,lat)+0.5_r8*divbz(k,i))*bz(i,lat)**2
981 : enddo ! k=kbot,nlev
982 : enddo ! i=i0,i1
983 168562680 : diag8(:,i0:i1,lat) = hdzmbz(:,i0:i1,lat) ! HDZMBZ
984 168562680 : diag9(:,i0:i1,lat) = hdzpbz(:,i0:i1,lat) ! HDZPBZ
985 : !
986 : ! Sum O+ at time n-1 to explicit terms: N(O+)/(2*DT) (N-1)
987 : ! 'EXPLICIT2' (zero at the poles)
988 : !
989 : !$omp parallel do private( i, k )
990 1393080 : do i=i0,i1
991 60545400 : do k=kbot,nlev
992 118304640 : explicit(k,i,lat) = explicit(k,i,lat)-(opnm_smooth(k,i,lat)-shapiro_const* &
993 118304640 : (opnm_smooth(k,i+2,lat)+opnm_smooth(k,i-2,lat)-4._r8* &
994 118304640 : (opnm_smooth(k,i+1,lat)+opnm_smooth(k,i-1,lat))+6._r8* &
995 354913920 : opnm_smooth(k,i,lat)))*dtx2inv
996 : op_dt(k,i,lat) = -(opnm_smooth(k,i,lat)-shapiro_const* &
997 : (opnm_smooth(k,i+2,lat)+opnm_smooth(k,i-2,lat)-4._r8* &
998 : (opnm_smooth(k,i+1,lat)+opnm_smooth(k,i-1,lat))+6._r8* &
999 60438240 : opnm_smooth(k,i,lat)))*dtx2inv
1000 : enddo ! k=kbot,nlev
1001 : enddo ! i=i0,i1
1002 168562680 : diag10(:,i0:i1,lat) = explicit(:,i0:i1,lat) ! EXPLICIT2
1003 : !
1004 : ! Begin coefficients p_coeff, q_coeff, r_coeff
1005 : !
1006 : !$omp parallel do private( i, k )
1007 1393080 : do i=i0,i1
1008 59259480 : do k=kbot,nlev-1
1009 57866400 : p_coeff(k,i,lat) = hdzmbz(k,i,lat)*djint(k ,i,lat)*tphdz0(k ,i,lat)
1010 57866400 : q_coeff(k,i,lat) = -(hdzpbz(k,i,lat)*djint(k+1,i,lat)*tphdz0(k+1,i,lat)+ &
1011 57866400 : hdzmbz(k,i,lat)*djint(k ,i,lat)*tphdz1(k ,i,lat))
1012 59152320 : r_coeff(k,i,lat) = hdzpbz(k,i,lat)*djint(k+1,i,lat)*tphdz1(k+1,i,lat)
1013 : enddo ! k=kbot,nlev-1
1014 : enddo ! i=i0,i1
1015 :
1016 168562680 : diag11(:,i0:i1,lat) = p_coeff(:,i0:i1,lat) ! P_COEFF0 (zero at poles)
1017 168562680 : diag12(:,i0:i1,lat) = q_coeff(:,i0:i1,lat) ! Q_COEFF0 (zero at ubc)
1018 168562680 : diag13(:,i0:i1,lat) = r_coeff(:,i0:i1,lat) ! R_COEFF0 (zero at ubc)
1019 : !
1020 : ! bdotu = B.U
1021 : ! Introducing neutral winds.
1022 : ! Am not using 0.5*(om(k)+om(k+1)) here because waccm omega is on midpoints (?)
1023 : ! (tiegcm has 0.5*(w(k,i,j0)+w(k+1,i,j0)) )
1024 : !
1025 : !$omp parallel do private( i, k )
1026 1393080 : do i=i0,i1
1027 60545400 : do k=kbot,nlev
1028 177456960 : bdotu(k,i,lat) = bx(i,lat)*un(k,i,lat)+by(i,lat)*vn(k,i,lat)+ &
1029 237895200 : hj(k,i,lat)*bz(i,lat)*om(k,i,lat)
1030 : enddo ! k=kbot,nlev
1031 : enddo ! i=i0,i1
1032 168562680 : diag14(:,i0:i1,lat) = bdotu(:,i0:i1,lat) ! BDOTU
1033 : !
1034 : ! Continue coefficients with vertical ion drift:
1035 : ! wi is converted from interfaces to midpoints (first use of wi).
1036 : ! The p,q,r coeffs are still zero at top boundary k=nlev, and at poles.
1037 : !
1038 : !$omp parallel do private( i, k )
1039 1393080 : do i=i0,i1
1040 57973560 : do k=kbot,nlev-2
1041 :
1042 169741440 : p_coeff(k+1,i,lat) = p_coeff(k+1,i,lat)+(bz(i,lat)*bdotu(k,i,lat)+ &
1043 226321920 : 0.5_r8*(wi(k+1,i,lat)+wi(k+2,i,lat)))*0.5_r8*hdz(k+1,i,lat)
1044 :
1045 56580480 : q_coeff(k,i,lat) = q_coeff(k,i,lat)-0.5_r8*(wi(k,i,lat)+wi(k+1,i,lat))*6._r8/Rearth
1046 :
1047 0 : r_coeff(k,i,lat) = r_coeff(k,i,lat)-(bz(i,lat)*bdotu(k+1,i,lat)+ &
1048 57866400 : 0.5_r8*(wi(k,i,lat)+wi(k+1,i,lat)))*0.5_r8*hdz(k,i,lat)
1049 :
1050 : enddo ! k=kbot,nlev-1
1051 : enddo ! i=i0,i1
1052 :
1053 168562680 : diag22(:,i0:i1,lat) = p_coeff(:,i0:i1,lat) ! P_COEFF0a
1054 168562680 : diag23(:,i0:i1,lat) = q_coeff(:,i0:i1,lat) ! Q_COEFF0a
1055 168562680 : diag24(:,i0:i1,lat) = r_coeff(:,i0:i1,lat) ! R_COEFF0a
1056 : !
1057 : ! Upper (nlev) and lower (kbot) boundaries of p,q,r_coeff:
1058 : ! (convert wi to midpoints)
1059 : !
1060 : ! tiegcm considers nlev-1 to be the top level. Do it tiegcm-style here,
1061 : ! and then extrapolate to nlev.
1062 : !
1063 : !$omp parallel do private( i )
1064 1393080 : do i=i0,i1
1065 2571840 : p_coeff(kbot,i,lat) = p_coeff(kbot,i,lat)+(bz(i,lat)* & ! reset p_coeff lbc
1066 1285920 : (2._r8*bdotu(kbot,i,lat)-bdotu(kbot+1,i,lat))+ &
1067 3857760 : 0.5_r8*(wi(kbot,i,lat)+wi(kbot+1,i,lat)))*0.5_r8*hdz(kbot,i,lat)
1068 :
1069 1285920 : q_coeff(nlev-1,i,lat) = q_coeff(nlev-1,i,lat)- &
1070 2571840 : 0.5_r8*(wi(nlev,i,lat)+wi(nlev-1,i,lat))*6._r8/Rearth
1071 :
1072 0 : r_coeff(nlev-1,i,lat) = r_coeff(nlev-1,i,lat)-(bz(i,lat)* &
1073 1285920 : (2._r8*bdotu(nlev-1,i,lat)-bdotu(nlev-2,i,lat))+ &
1074 2679000 : 0.5_r8*(wi(nlev,i,lat)+wi(nlev-1,i,lat)))*0.5_r8*hdz(nlev-1,i,lat)
1075 : enddo ! i=i0,i1
1076 : !
1077 : ! Extrapolate to top level (tiegcm does not do this):
1078 : !
1079 214320 : p_coeff(nlev,i0:i1,lat) = 1.5_r8*p_coeff(nlev-1,i0:i1,lat)- &
1080 1607400 : 0.5_r8*p_coeff(nlev-2,i0:i1,lat)
1081 : q_coeff(nlev,i0:i1,lat) = 1.5_r8*q_coeff(nlev-1,i0:i1,lat)- &
1082 1393080 : 0.5_r8*q_coeff(nlev-2,i0:i1,lat)
1083 : r_coeff(nlev,i0:i1,lat) = 1.5_r8*r_coeff(nlev-1,i0:i1,lat)- &
1084 1393080 : 0.5_r8*r_coeff(nlev-2,i0:i1,lat)
1085 : !
1086 : ! All P,Q,R are zero at the poles. Polar values will be set after third lat scan.
1087 168562680 : diag15(:,i0:i1,lat) = p_coeff(:,i0:i1,lat) ! P_COEFF1 (zero at ubc and poles)
1088 168562680 : diag17(:,i0:i1,lat) = r_coeff(:,i0:i1,lat) ! R_COEFF1 (ok at ubc, zero at poles)
1089 : !
1090 : ! Additions to Q coefficients (includes q_coeff lbc,ubc):
1091 : !$omp parallel do private( i, k )
1092 1393080 : do i=i0,i1
1093 60545400 : do k=kbot,nlev
1094 60438240 : q_coeff(k,i,lat) = q_coeff(k,i,lat)-bdotu(k,i,lat)*dvb(i,lat)*bz(i,lat)-dtx2inv
1095 : enddo ! k=kbot,nlev-1
1096 : enddo ! i=i0,i1
1097 : !
1098 : ! Plot Q_COEFF1 after ubc has been set.
1099 168562680 : diag16(:,i0:i1,lat) = q_coeff(:,i0:i1,lat) ! Q_COEFF1 (ok at ubc, zero at poles)
1100 : !
1101 : ! Upper boundary condition for O+:
1102 : !$omp parallel do private( i )
1103 1393080 : do i=i0,i1
1104 1285920 : ubca(i) = 0._r8
1105 1285920 : ubcb(i) = -bz(i,lat)**2*djint(nlev,i,lat)*tphdz0(nlev,i,lat)-ubca(i)
1106 1285920 : ubca(i) = -bz(i,lat)**2*djint(nlev,i,lat)*tphdz1(nlev,i,lat)+ubca(i)
1107 : !
1108 : ! Q = Q+B/A*R
1109 : q_coeff(nlev,i,lat) = q_coeff(nlev,i,lat)+ubcb(i)/ubca(i)* &
1110 1285920 : r_coeff(nlev,i,lat)
1111 : !
1112 : ! F = F -R/A*PHI
1113 : explicit(nlev,i,lat) = explicit(nlev,i,lat)-opflux(i,lat)* & ! explicit ubc
1114 1285920 : r_coeff(nlev,i,lat)/ubca(i)
1115 1393080 : r_coeff(nlev,i,lat) = 0._r8 ! r_coeff ubc is reset to zero
1116 : enddo ! i=i0,i1
1117 : !
1118 : ! Ubc of EXPLICIT3 has a stripe along the mag equator, unlike the level below.
1119 : !
1120 168562680 : diag18(:,i0:i1,lat) = explicit(:,i0:i1,lat) ! EXPLICIT3 (ubc ok, zero at poles)
1121 168562680 : diag19(:,i0:i1,lat) = p_coeff(:,i0:i1,lat) ! P_COEFF2 (zero at ubc, zero at poles)
1122 168562680 : diag20(:,i0:i1,lat) = q_coeff(:,i0:i1,lat) ! Q_COEFF2 (ubc ok, zero at poles)
1123 168562680 : diag21(:,i0:i1,lat) = r_coeff(:,i0:i1,lat) ! R_COEFF2 (zero at ubc, zero at poles)
1124 : !
1125 : ! At this point, TIEGCM calculates "sources and sinks" xiop2p and xiop2d.
1126 : ! Also calculates op_loss, which is subtracted from q_coeff.
1127 : ! Then TIEGCM "Add source term to RHS (explicit terms)", and calculates
1128 : ! lower boundary condition N(O+) = Q/L (q_coeff, explicit, p_coeff), and
1129 : ! finally calls trsolv.
1130 : !
1131 36480 : 300 continue
1132 : enddo ! end third latitude scan (lat=lat0,lat1)
1133 : !
1134 : !------------------------ End third latitude scan ---------------------
1135 :
1136 : !
1137 : ! Set poles for selected diagnostics:
1138 : !
1139 123703680 : call setpoles(diag26(kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1) ! EXPLICITa
1140 123703680 : call setpoles(diag27(kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1) ! EXPLICITb
1141 123703680 : call setpoles(diag2 (kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1) ! EXPLICIT0
1142 123703680 : call setpoles(diag3 (kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1) ! EXPLICIT1
1143 123703680 : call setpoles(diag6 (kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1) ! DJINT
1144 : !
1145 : ! All tasks have global 2d bmod2.
1146 : ! bmod2 was set by sub magfield (getapex.F90)
1147 : ! allocate(bmod2(0:nlonp1,jspole-1:jnpole+1))
1148 : ! Copy bmod2 poles to diagnostic array.
1149 : !
1150 : !$omp parallel do private( i, k )
1151 474240 : do i=i0,i1
1152 20611200 : do k=kbot,nlev
1153 20136960 : diag25(k,i,j0) = bmod2(i,j0)
1154 20574720 : diag25(k,i,j1) = bmod2(i,j1)
1155 : enddo
1156 : enddo
1157 36480 : if (debug_hist) then
1158 0 : call savefld_waccm(diag25,'BMOD2' ,nlev,i0,i1,j0,j1)
1159 : endif
1160 : !
1161 : ! Assign polar values to coefficients for trsolv.
1162 : !
1163 123703680 : call setpoles(explicit(kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
1164 123703680 : call setpoles(p_coeff (kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
1165 123703680 : call setpoles(q_coeff (kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
1166 123703680 : call setpoles(r_coeff (kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
1167 : !
1168 : ! Call solver, defining O+ output op_out:
1169 : !
1170 : ! Its best not to call this unless the coefficients and explicit terms
1171 : ! have been properly set in the third latitude scan above (e.g., during
1172 : ! "goto 300" debugging above, where the coeffs may not have been calculated).
1173 : !
1174 : if (calltrsolv) then
1175 :
1176 : !$omp parallel do private( lat )
1177 145920 : do lat=j0,j1
1178 :
1179 109440 : call trsolv(p_coeff (kbot:nlev,i0:i1,lat), &
1180 109440 : q_coeff (kbot:nlev,i0:i1,lat), &
1181 109440 : r_coeff (kbot:nlev,i0:i1,lat), &
1182 109440 : explicit(kbot:nlev,i0:i1,lat), &
1183 109440 : op_out (kbot:nlev,i0:i1,lat), &
1184 309277440 : kbot,nlev,kbot,nlev,i0,i1 )
1185 :
1186 :
1187 :
1188 :
1189 : ! for term analysis
1190 5070720 : do k=kbot,nlev-1
1191 : ! diffusion
1192 9849600 : amb_diff(k,i0:i1,lat) = amb_diff(k,i0:i1,lat) + &
1193 4924800 : hdzmbz(k,i0:i1,lat)*djint(k, i0:i1,lat)*tphdz0(k, i0:i1,lat)* op_out(k-1,i0:i1,lat) &
1194 4924800 : -(hdzpbz(k,i0:i1,lat)*djint(k+1,i0:i1,lat)*tphdz0(k+1,i0:i1,lat)+ &
1195 : hdzmbz(k,i0:i1,lat)*djint(k ,i0:i1,lat)*tphdz1(k ,i0:i1,lat))* op_out(k,i0:i1,lat) &
1196 83721600 : +hdzpbz(k,i0:i1,lat)*djint(k+1,i0:i1,lat)*tphdz1(k+1,i0:i1,lat)* op_out(k+1,i0:i1,lat)
1197 :
1198 : ! electric field transport
1199 4924800 : if (k <= nlev-2) then
1200 : dfield(k,i0:i1,lat) = dfield(k,i0:i1,lat)+ &
1201 4815360 : (0.5_r8*(wi(k+1,i0:i1,lat)+wi(k+2,i0:i1,lat)))* &
1202 : 0.5_r8*hdz(k+1,i0:i1,lat)* op_out(k-1,i0:i1,lat) &
1203 : -0.5_r8*(wi(k,i0:i1,lat)+wi(k+1,i0:i1,lat))* &
1204 : 6._r8/Rearth* op_out(k,i0:i1,lat) &
1205 : -(0.5_r8*(wi(k,i0:i1,lat)+wi(k+1,i0:i1,lat)))*0.5_r8*hdz(k,i0:i1,lat) &
1206 67415040 : * op_out(k+1,i0:i1,lat)
1207 : else
1208 : dfield(k,i0:i1,lat) = dfield(k,i0:i1,lat)+ &
1209 : (1*(wi(k+1,i0:i1,lat)))* &
1210 : 0.5_r8*hdz(k+1,i0:i1,lat)* op_out(k-1,i0:i1,lat) &
1211 : -0.5_r8*(wi(k,i0:i1,lat)+wi(k+1,i0:i1,lat))* &
1212 : 6._r8/Rearth* op_out(k,i0:i1,lat) &
1213 : -(0.5_r8*(wi(k,i0:i1,lat)+wi(k+1,i0:i1,lat)))*0.5_r8*hdz(k,i0:i1,lat) &
1214 1422720 : * op_out(k+1,i0:i1,lat)
1215 : endif
1216 : ! wind transport
1217 : dwind(k,i0:i1,lat)= &
1218 0 : (bz(i0:i1,lat)*bdotu(k,i0:i1,lat))* 0.5_r8*hdz(k+1,i0:i1,lat) * op_out(k-1,i0:i1,lat) &
1219 : -bdotu(k,i0:i1,lat)*dvb(i0:i1,lat)*bz(i0:i1,lat)* op_out(k,i0:i1,lat) &
1220 64022400 : -(bz(i0:i1,lat)*bdotu(k+1,i0:i1,lat))*0.5_r8*hdz(k,i0:i1,lat)* op_out(k+1,i0:i1,lat)
1221 :
1222 : ! dO+/dt
1223 64131840 : op_dt(k,i0:i1,lat)= dtx2inv* op_out(k,i0:i1,lat) + op_dt(k,i0:i1,lat)
1224 : !
1225 : enddo ! k=lev0+1,lev1-1
1226 :
1227 : enddo
1228 :
1229 36480 : call savefld_waccm(op_dt,'op_dt',nlev,i0,i1,j0,j1)
1230 36480 : call savefld_waccm(amb_diff,'amb_diff',nlev,i0,i1,j0,j1)
1231 36480 : call savefld_waccm(dfield,'dfield',nlev,i0,i1,j0,j1)
1232 36480 : call savefld_waccm(dwind,'dwind',nlev,i0,i1,j0,j1)
1233 :
1234 36480 : if (debug_hist) then
1235 0 : call savefld_waccm(op_out,'OP_SOLVE',nlev,i0,i1,j0,j1)
1236 : endif
1237 :
1238 : else ! trsolv not called (debug only)
1239 : op_out (kbot:nlev,i0:i1,j0:j1) = op (kbot:nlev,i0:i1,j0:j1)
1240 : opnm_out(kbot:nlev,i0:i1,j0:j1) = opnm(kbot:nlev,i0:i1,j0:j1)
1241 : endif ! calltrsolv
1242 : !
1243 : ! Write fields from third latitude scan to waccm history:
1244 : !
1245 36480 : if (debug_hist) then
1246 0 : call savefld_waccm(explicit,'EXPLICIT',nlev,i0,i1,j0,j1) ! non-zero at ubc
1247 0 : call savefld_waccm(p_coeff ,'P_COEFF' ,nlev,i0,i1,j0,j1) ! zero at ubc?
1248 0 : call savefld_waccm(q_coeff ,'Q_COEFF' ,nlev,i0,i1,j0,j1) ! non-zero at ubc
1249 0 : call savefld_waccm(r_coeff ,'R_COEFF' ,nlev,i0,i1,j0,j1) ! is set zero at ubc
1250 :
1251 0 : call savefld_waccm(bdotdh_diff(:,i0:i1,j0:j1), 'BDOTDH_DIFF',nlev,i0,i1,j0,j1)
1252 0 : call savefld_waccm(diag1 ,'BDZDVB_OP',nlev,i0,i1,j0,j1)
1253 0 : call savefld_waccm(diag2 ,'EXPLICIT0',nlev,i0,i1,j0,j1)
1254 0 : call savefld_waccm(diag26,'EXPLICITa',nlev,i0,i1,j0,j1)
1255 0 : call savefld_waccm(diag27,'EXPLICITb',nlev,i0,i1,j0,j1)
1256 0 : call savefld_waccm(diag3 ,'EXPLICIT1',nlev,i0,i1,j0,j1)
1257 0 : call savefld_waccm(diag4 ,'TPHDZ0' ,nlev,i0,i1,j0,j1)
1258 0 : call savefld_waccm(diag5 ,'TPHDZ1' ,nlev,i0,i1,j0,j1)
1259 0 : call savefld_waccm(diag6 ,'DJINT' ,nlev,i0,i1,j0,j1)
1260 0 : call savefld_waccm(diag7 ,'DIVBZ' ,nlev,i0,i1,j0,j1)
1261 0 : call savefld_waccm(diag8 ,'HDZMBZ' ,nlev,i0,i1,j0,j1)
1262 0 : call savefld_waccm(diag9 ,'HDZPBZ' ,nlev,i0,i1,j0,j1)
1263 0 : call savefld_waccm(diag10,'EXPLICIT2',nlev,i0,i1,j0,j1)
1264 0 : call savefld_waccm(diag11,'P_COEFF0' ,nlev,i0,i1,j0,j1)
1265 0 : call savefld_waccm(diag12,'Q_COEFF0' ,nlev,i0,i1,j0,j1)
1266 0 : call savefld_waccm(diag13,'R_COEFF0' ,nlev,i0,i1,j0,j1)
1267 0 : call savefld_waccm(diag14,'BDOTU' ,nlev,i0,i1,j0,j1)
1268 0 : call savefld_waccm(diag15,'P_COEFF1' ,nlev,i0,i1,j0,j1)
1269 0 : call savefld_waccm(diag16,'Q_COEFF1' ,nlev,i0,i1,j0,j1)
1270 0 : call savefld_waccm(diag17,'R_COEFF1' ,nlev,i0,i1,j0,j1)
1271 0 : call savefld_waccm(diag18,'EXPLICIT3',nlev,i0,i1,j0,j1)
1272 0 : call savefld_waccm(diag19,'P_COEFF2' ,nlev,i0,i1,j0,j1)
1273 0 : call savefld_waccm(diag20,'Q_COEFF2' ,nlev,i0,i1,j0,j1)
1274 0 : call savefld_waccm(diag21,'R_COEFF2' ,nlev,i0,i1,j0,j1)
1275 0 : call savefld_waccm(diag22,'P_COEFF0a',nlev,i0,i1,j0,j1)
1276 0 : call savefld_waccm(diag23,'Q_COEFF0a',nlev,i0,i1,j0,j1)
1277 0 : call savefld_waccm(diag24,'R_COEFF0a',nlev,i0,i1,j0,j1)
1278 : endif
1279 : !
1280 : !------------------------------------------------------------------------
1281 : !
1282 : ! Filter O+ output from solver:
1283 : ! (TIMEGCM calls both filters, whereas TIEGCM calls only filter2)
1284 : !
1285 : ! call filter1_op(op_out(kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
1286 : !
1287 36480 : if (ring_polar_filter) then
1288 123703680 : call ringfilter_op(op_out(kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
1289 : else
1290 0 : call filter2_op(op_out(kbot:nlev,i0:i1,j0:j1),kbot,nlev,i0,i1,j0,j1)
1291 : endif
1292 : !
1293 : !----------------------- Begin fourth latitude scan ---------------------
1294 : !
1295 : !$omp parallel do private(lat, i, k, opfloor)
1296 145920 : do lat=j0,j1
1297 1422720 : do i=i0,i1
1298 61833600 : do k=kbot,nlev
1299 181232640 : opnm_out(k,i,lat) = dtsmooth*op(k,i,lat)+dtsmooth_div2* &
1300 242956800 : (opnm(k,i,lat)+op_out(k,i,lat))
1301 : enddo
1302 : enddo
1303 : !
1304 : ! Insure non-negative O+ output:
1305 1422720 : do i=i0,i1
1306 61833600 : do k=kbot,nlev
1307 60410880 : if (op_out (k,i,lat) < 1.e-5_r8) op_out (k,i,lat) = 1.e-5_r8
1308 61724160 : if (opnm_out(k,i,lat) < 1.e-5_r8) opnm_out(k,i,lat) = 1.e-5_r8
1309 : enddo ! k=lev0,lev1-1
1310 : enddo ! i=lon0,lon1
1311 : !
1312 : ! Enforce O+ minimum if enforce_opfloor is true.
1313 : ! Opfloor is Stan's "smooth floor" (product of two Gaussians,
1314 : ! dependent on latitude and pressure level) (opmin=3000.0):
1315 : !
1316 145920 : if (enforce_floor) then
1317 5143680 : zpmid(kbot:nlev) = log(50.e-6_r8/pmid(kbot:nlev)) ! tgcm levs -- maybe done once at init time
1318 5143680 : do k=kbot,nlev
1319 0 : opfloor = opmin*exp(-(glat(lat)/90.0_r8)**2/0.3_r8) &
1320 5034240 : *exp(-((zpmid(k)-4.25_r8)/zpmid(nlev))**2/0.1_r8)
1321 65554560 : do i=i0,i1
1322 65445120 : if (op_out(k,i,lat) < opfloor) then
1323 1412235 : op_out(k,i,lat) = opfloor
1324 : endif ! opout < opfloor
1325 : enddo ! i=lon0,lon1
1326 : enddo ! k=lev0,lev1-1
1327 : endif ! enforce_opfloor
1328 :
1329 : enddo ! lat=lat0,lat1
1330 :
1331 : !
1332 : ! Save O+ output to WACCM history (cm^3):
1333 36480 : if (debug_hist) then
1334 0 : call savefld_waccm(op_out (:,i0:i1,j0:j1),'OP_OUT' ,nlev,i0,i1,j0,j1)
1335 0 : call savefld_waccm(opnm_out(:,i0:i1,j0:j1),'OPNM_OUT',nlev,i0,i1,j0,j1)
1336 : endif
1337 36480 : end subroutine oplus_xport
1338 : !-----------------------------------------------------------------------
1339 36480 : subroutine oplus_flux(opflux,lon0,lon1,lat0,lat1)
1340 : !
1341 : ! Calculate O+ number flux for sub oplus_xport.
1342 : ! Flux is returned in opflux(lon0:lon1,lat0:lat1).
1343 : !
1344 : ! alatm: geomagnetic latitude at each geographic grid point (radians)
1345 : use getapex,only: alatm ! (nlonp1,jspole:jnpole)
1346 : !
1347 : ! Args:
1348 : integer,intent(in) :: lon0,lon1,lat0,lat1
1349 : real(r8),intent(out) :: opflux(lon0:lon1,lat0:lat1)
1350 : !
1351 : ! Local:
1352 : integer :: i,j
1353 72960 : real(r8),dimension(lon0:lon1,lat0:lat1) :: chi ! solar zenith angle
1354 : real(r8),parameter :: &
1355 : phid = 2.0e8_r8, &
1356 : phin = -2.0e8_r8, &
1357 : ! phin = 0._r8, &
1358 : ppolar = 0._r8
1359 72960 : real(r8) :: a(lon0:lon1)
1360 72960 : real(r8) :: fed(lon0:lon1)
1361 72960 : real(r8) :: fen(lon0:lon1)
1362 : !
1363 : ! Set some paramaters:
1364 36480 : pi = 4._r8*atan(1._r8)
1365 36480 : rtd = 45._r8/atan(1._r8)
1366 : !
1367 : ! Sub get_zenith calls sub zenith (..cam/src/physics/cam/zenith.F90)
1368 36480 : call get_zenith(chi,lon0,lon1,lat0,lat1)
1369 : !
1370 : ! Latitude scan:
1371 145920 : do j=lat0,lat1
1372 : !
1373 : ! Longitude loop:
1374 1459200 : do i=lon0,lon1
1375 1313280 : if (abs(alatm(i,j))-pi/24._r8>=0._r8) then
1376 1208495 : a(i) = 1._r8
1377 : else
1378 104785 : a(i)=.5_r8*(1._r8+sin(pi*(abs(alatm(i,j))-pi/48._r8)/(pi/24._r8)))
1379 104785 : if (a(i) < 0.05_r8) a(i) = 0.05_r8
1380 : endif
1381 1313280 : fed(i) = phid*a(i)
1382 1313280 : fen(i) = phin*a(i)
1383 1313280 : if (chi(i,j)-0.5_r8*pi >= 0._r8) then
1384 656640 : opflux(i,j) = fen(i)
1385 : else
1386 656640 : opflux(i,j) = fed(i)
1387 : endif
1388 1313280 : if ((chi(i,j)*rtd-80._r8)*(chi(i,j)*rtd-100._r8) < 0._r8) then
1389 : opflux(i,j) = .5_r8*(fed(i)+fen(i))+.5_r8*(fed(i)-fen(i))* &
1390 220860 : cos(pi*(chi(i,j)*rtd-80._r8)/20._r8)
1391 : endif
1392 : !
1393 : ! Add ppolar if magnetic latitude >= 60 degrees:
1394 : ! QUESTION: is the 60 deg here related to critical angles crit(2) in tiegcm?
1395 : ! 3/15/15: opflux is comparable to tiegcm.
1396 : !
1397 1313280 : if (abs(alatm(i,j))-pi/3._r8 >= 0._r8) &
1398 540265 : opflux(i,j) = opflux(i,j)+ppolar
1399 : enddo ! i=lon0,lon1
1400 : enddo ! j=lat0,lat1
1401 : !
1402 36480 : end subroutine oplus_flux
1403 : !-----------------------------------------------------------------------
1404 72960 : subroutine get_zenith(chi,i0,i1,j0,j1)
1405 : !
1406 : ! Get solar zenith angle chi(i0:i1,j0:j1) (radians)
1407 : ! Subroutine zenith returns cos(chi) at each (i,j)
1408 : ! Note glon(i0:i1) from edyn_init is in -180 -> +180 (TIEGCM mode)
1409 : !
1410 : use time_manager,only : get_curr_calday
1411 : use edyn_geogrid,only : glon,glat
1412 : use orbit, only : zenith
1413 : !
1414 : ! Args:
1415 : integer,intent(in) :: i0,i1,j0,j1
1416 : real(r8),intent(out) :: chi(i0:i1,j0:j1)
1417 : !
1418 : ! Local:
1419 : integer :: i,j
1420 : real(r8) :: dtr,calday
1421 : real(r8) :: cosZenAngR(1)
1422 :
1423 36480 : dtr = pi/180._r8
1424 36480 : calday = get_curr_calday() ! fractional day of year
1425 145920 : do j=j0,j1
1426 1459200 : do i=i0,i1
1427 3939840 : call zenith(calday,(/dtr*glat(j)/),(/dtr*glon(i)/),cosZenAngR,1)
1428 1422720 : chi(i,j) = acos(cosZenAngR(1))
1429 : enddo
1430 : enddo
1431 36480 : end subroutine get_zenith
1432 : !-----------------------------------------------------------------------
1433 36480 : subroutine divb(dvb,i0,i1,j0,j1)
1434 : !
1435 : ! Evaluate divergence of B, the unit magnetic field vector.
1436 : ! (all processors have the full global 2d field)
1437 : !
1438 : ! Args:
1439 : integer,intent(in) :: i0,i1,j0,j1
1440 : real(r8),intent(out) :: dvb(i0:i1,j0:j1)
1441 : !
1442 : ! Local:
1443 : integer :: i,j,jm1,jp1
1444 :
1445 1459200 : dvb = 0._r8
1446 :
1447 36480 : if (debug_hist) then
1448 0 : call savefld_waccm(bx(i0:i1,j0:j1),'OPLUS_BX',1,i0,i1,j0,j1)
1449 0 : call savefld_waccm(by(i0:i1,j0:j1),'OPLUS_BY',1,i0,i1,j0,j1)
1450 0 : call savefld_waccm(bz(i0:i1,j0:j1),'OPLUS_BZ',1,i0,i1,j0,j1)
1451 0 : call savefld_waccm(bmod2(i0:i1,j0:j1),'OPLUS_BMAG',1,i0,i1,j0,j1)
1452 : endif
1453 : !
1454 : ! Note Rearth is in cm.
1455 : ! (bx,by,bz are set by sub magfield (getapex.F90))
1456 : ! (dphi,dlamda, and cs are set by sub set_geogrid (edyn_init.F90))
1457 : !
1458 145920 : do j=j0,j1
1459 109440 : jm1 = j-1
1460 109440 : jp1 = j+1
1461 1459200 : do i=i0,i1
1462 3939840 : dvb(i,j) = (((bx(i+1,j)-bx(i-1,j))/(2._r8*dlamda)+ &
1463 1313280 : (cs(jp1)*by(i,jp1)-cs(jm1)*by(i,jm1))/(2._r8*dphi))/ &
1464 6675840 : cs(j)+2._r8*bz(i,j))/Rearth
1465 : enddo ! i=i0,i1
1466 : enddo ! j=j0,j1
1467 36480 : end subroutine divb
1468 : !-----------------------------------------------------------------------
1469 321480 : subroutine rrk(t,rms,ps1,ps2,n2,tr,ans,lon0,lon1,lev0,lev1)
1470 : !
1471 : ! Returns ambipolar diffusion coefficient in ans.
1472 : !
1473 : ! Args:
1474 : integer,intent(in) :: lon0,lon1,lev0,lev1
1475 : real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: &
1476 : t,rms,ps1,ps2,n2,tr
1477 : real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: ans
1478 : !
1479 : ! Local:
1480 : !
1481 : integer :: k,i
1482 : !
1483 : !$omp parallel do private(i,k)
1484 4179240 : do i=lon0,lon1
1485 177456960 : do k=lev0,lev1-1
1486 :
1487 347198400 : ans(k,i) = 1.42e17_r8*boltz*t(k,i)/(p0*expz(k)*.5_r8*(rms(k,i)+ &
1488 173599200 : rms(k+1,i))*(ps2(k,i)*rmassinv_o1*sqrt(tr(k,i))*(1._r8-0.064_r8* &
1489 : log10(tr(k,i)))**2*colfac+18.6_r8*n2(k,i)*rmassinv_n2+18.1_r8* &
1490 698254560 : ps1(k,i)*rmassinv_o2))
1491 :
1492 : enddo ! k=lev0,lev1
1493 4179240 : ans(lev1,i) = ans(lev1-1,i) ! should not need to do this
1494 :
1495 : enddo ! i=lon0,lon1
1496 : !
1497 : ! Cap ambipolar diffusion coefficient in ans.
1498 : !
1499 : ! acceptable range for limiter 1.e8 to 1.e9 ...
1500 181636200 : where( ans(:,:) > adiff_limiter )
1501 : ans(:,:) = adiff_limiter
1502 : endwhere
1503 :
1504 321480 : end subroutine rrk
1505 : !-----------------------------------------------------------------------
1506 321480 : subroutine diffus(tp,en,hj,ans,i0,i1,lev0,lev1,lat)
1507 : ! kbot,nlev
1508 : ! Evaluates ans = (d/(h*dz)*tp+m*g/r)*en
1509 : ! Remember: "bot2top": lev0=kbot=bottom, lev1=nlev=top
1510 : !
1511 : ! Args:
1512 : integer :: i0,i1,lev0,lev1,lat
1513 : real(r8),dimension(lev0:lev1,i0:i1),intent(in) :: tp,en,hj
1514 : real(r8),dimension(lev0:lev1,i0:i1),intent(out) :: ans
1515 : !
1516 : ! Local:
1517 : integer :: i,k
1518 : real(r8) :: mgr
1519 :
1520 321480 : mgr = rmass_op*grav_cm/gask
1521 :
1522 : !$omp parallel do private(i,k)
1523 4179240 : do i=i0,i1
1524 173599200 : do k=lev0,lev1-2
1525 509224320 : ans(k+1,i) = 1._r8/(2._r8*hj(k+1,i)*dzp)*(tp(k+2,i)*en(k+2,i)- &
1526 682823520 : tp(k,i)*en(k,i))+mgr*en(k+1,i)
1527 : enddo
1528 321480 : if (debug) then
1529 : write(iulog,"('diffus: lat=',i4,' i=',i4,' ans(lev0:lev1-1,i)=',2es12.4)") &
1530 : lat,i,minval(ans(lev0:lev1-1,i)),maxval(ans(lev0:lev1-1,i))
1531 : endif
1532 : enddo
1533 : !
1534 : ! Upper and lower boundaries:
1535 : !
1536 : !$omp parallel do private(i)
1537 4179240 : do i=i0,i1
1538 : !
1539 : ! Upper boundary:
1540 7715520 : ans(lev1,i) = 1._r8/(hj(lev1,i)*dzp)*(tp(lev1,i)*en(lev1,i)- &
1541 11573280 : tp(lev1-1,i)*en(lev1-1,i))+mgr*en(lev1,i)
1542 : !
1543 : ! Lower boundary:
1544 7715520 : ans(lev0,i) = 1._r8/(hj(lev0,i)*dzp)*(tp(lev0+1,i)*en(lev0+1,i)- &
1545 11894760 : tp(lev0,i)*en(lev0,i))+mgr*en(lev0,i)
1546 : enddo
1547 321480 : end subroutine diffus
1548 : !-----------------------------------------------------------------------
1549 535800 : subroutine bdotdh(phijm1,phij,phijp1,ans,lon0,lon1,lev0,lev1,lat)
1550 : !
1551 : ! Evaluates ans = (b(h)*del(h))*phi
1552 : !
1553 : ! Args:
1554 : integer,intent(in) :: lon0,lon1,lev0,lev1,lat
1555 : real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: phijm1,phijp1
1556 : real(r8),dimension(lev0:lev1,lon0-2:lon1+2),intent(inout) :: phij ! why intent(inout)?
1557 : real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: ans
1558 : !
1559 : ! Local:
1560 : integer :: k,i
1561 : !
1562 : ! Note phij longitude dimension is lon0-2:lon1+2 (only i-1 and i+1 are used).
1563 : ! Halo longitudes i-1 and i+1 must have been set before this routine is
1564 : ! called. ('by' is use-associated above)
1565 : !
1566 : !$omp parallel do private( i, k )
1567 6965400 : do i=lon0,lon1
1568 302727000 : do k=lev0,lev1
1569 591523200 : ans(k,i) = 1._r8/Rearth*(bx(i,lat)/(cs(lat)*2._r8*dlamda)* &
1570 591523200 : (phij(k,i+1)-phij(k,i-1))+by(i,lat)* &
1571 1485237600 : (phijp1(k,i)-phijm1(k,i))/(2._r8*dphi))
1572 : enddo ! k=lev0,lev1
1573 : enddo ! i=lon0,lon1
1574 : !
1575 535800 : end subroutine bdotdh
1576 : !-----------------------------------------------------------------------
1577 107160 : subroutine bdzdvb(phi,dvb,h,ans,lev0,lev1,lon0,lon1,lat)
1578 : !
1579 : ! Evaluates ans = (bz*d/(h*dz)+divb)*phi
1580 : !
1581 : ! Args:
1582 : integer,intent(in) :: lev0,lev1,lon0,lon1,lat
1583 : real(r8),intent(in) :: dvb(lon0:lon1)
1584 : real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: phi,h
1585 : real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: ans
1586 : !
1587 : ! Local:
1588 : integer :: k,i
1589 : !
1590 : !$omp parallel do private( i, k )
1591 1393080 : do i=lon0,lon1
1592 57973560 : do k=lev0+1,lev1-1
1593 226321920 : ans(k,i) = bz(i,lat)/(2._r8*h(k,i)*dzp)*(phi(k+1,i)-phi(k-1,i))+ &
1594 284188320 : dvb(i)*phi(k,i)
1595 : enddo ! k=lev0+1,lev1-1
1596 : enddo ! i=lon0,lon1
1597 : !
1598 : ! Upper and lower boundaries:
1599 : !$omp parallel do private( i )
1600 1393080 : do i=lon0,lon1
1601 1285920 : ans(lev1,i) = bz(i,lat)/(h(lev1,i)*dzp)*(phi(lev1,i)- &
1602 2571840 : phi(lev1-1,i))+dvb(i)*phi(lev1,i)
1603 0 : ans(lev0,i) = bz(i,lat)/(h(lev0,i)*dzp)* &
1604 1393080 : (phi(lev0+1,i)-phi(lev0,i))+dvb(i)*phi(lev0,i)
1605 : enddo ! i=lon0,lon1
1606 107160 : end subroutine bdzdvb
1607 : !-----------------------------------------------------------------------
1608 : subroutine printpoles(f,klev,k0,k1,i0,i1,j0,j1,name)
1609 : use edyn_geogrid,only : nlat
1610 : !
1611 : ! Args:
1612 : integer,intent(in) :: klev,k0,k1,i0,i1,j0,j1
1613 : real(r8),intent(in) :: f(k0:k1,i0:i1,j0:j1)
1614 : character(len=*),intent(in) :: name
1615 : !
1616 : ! Print values at the poles at klev:
1617 : if (j0==1) then
1618 : if (debug.and.masterproc) write(iulog,"(/,'printpoles ',a,' spole: klev=',i4,' f(klev,i0:i1,j0)=',/,(8es12.4))") &
1619 : name,klev,f(klev,i0:i1,j0)
1620 : endif
1621 : if (j1==nlat) then
1622 : if (debug.and.masterproc) write(iulog,"(/,'printpoles ',a,' npole: klev=',i4,' f(klev,i0:i1,j1)=',/,(8es12.4))") &
1623 : name,klev,f(klev,i0:i1,j1)
1624 : endif
1625 :
1626 : end subroutine printpoles
1627 : !-----------------------------------------------------------------------
1628 : subroutine filter1_op(f,k0,k1,i0,i1,j0,j1)
1629 : !
1630 : ! Polar fft filter, option 1 (see filter.F90).
1631 : !
1632 : use filter_module,only: filter1,kut1
1633 : use edyn_mpi ,only: mp_gatherlons_f3d,mytidi
1634 : use edyn_mpi ,only: mp_scatterlons_f3d
1635 : use edyn_geogrid ,only: nlon
1636 : !
1637 : ! Args:
1638 : integer,intent(in) :: k0,k1,i0,i1,j0,j1
1639 : real(r8),intent(inout) :: f(k0:k1,i0:i1,j0:j1)
1640 : !
1641 : ! Local:
1642 : integer :: i,j,k,nlevs
1643 : real(r8) :: fik(nlon,k1-k0+1)
1644 : type(array_ptr_type) :: fkij(1) ! fkij(1)%ptr(k1-k0+1,nlon,j0:j1)
1645 :
1646 : nlevs = k1-k0+1
1647 : !
1648 : ! Define lons in fkij from current task subdomain:
1649 : !
1650 : allocate(fkij(1)%ptr(nlevs,nlon,j0:j1))
1651 : do j=j0,j1
1652 : do i=i0,i1
1653 : do k=k0,k1 ! kbot,nlev
1654 : fkij(1)%ptr(k-k0+1,i,j) = f(k,i,j)
1655 : enddo
1656 : enddo
1657 : enddo
1658 : !
1659 : ! Gather longitudes into tasks in first longitude column of task table
1660 : ! (leftmost of each j-row) for global fft. (i.e., tasks with mytidi==0
1661 : ! gather lons from other tasks in that row). This includes all latitudes.
1662 : !
1663 : call mp_gatherlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1)
1664 : !
1665 : ! Only leftmost tasks at each j-row of tasks does the global filtering:
1666 : !
1667 : if (mytidi==0) then
1668 : !
1669 : ! Define 2d array with all longitudes for filter at each latitude:
1670 : !
1671 : latscan: do j=j0,j1
1672 : if (kut1(j) >= nlon/2) cycle latscan
1673 : do i=1,nlon
1674 : do k=k0,k1
1675 : fik(i,k-k0+1) = fkij(1)%ptr(k-k0+1,i,j)
1676 : enddo
1677 : enddo
1678 : !
1679 : ! Remove wave numbers > kut(lat):
1680 : !
1681 : call filter1(fik,1,nlevs,j)
1682 : !
1683 : ! Return filtered array to fkij:
1684 : !
1685 : do i=1,nlon
1686 : do k=k0,k1
1687 : fkij(1)%ptr(k-k0+1,i,j) = fik(i,k-k0+1)
1688 : enddo
1689 : enddo ! i=1,nlon
1690 : enddo latscan ! j=j0,j1
1691 : endif ! mytidi==0
1692 : !
1693 : ! Now leftmost task at each j-row must redistribute filtered data
1694 : ! back to other tasks in the j-row (mytidi>0,mytidj) (includes latitude):
1695 : !
1696 : call mp_scatterlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1)
1697 : !
1698 : ! Return filtered array to inout field at task subdomain:
1699 : !
1700 : do j=j0,j1
1701 : do i=i0,i1
1702 : do k=k0,k1
1703 : f(k,i,j) = fkij(1)%ptr(k-k0+1,i,j)
1704 : enddo
1705 : enddo
1706 : enddo
1707 : deallocate(fkij(1)%ptr)
1708 : end subroutine filter1_op
1709 : !-----------------------------------------------------------------------
1710 0 : subroutine filter2_op(f,k0,k1,i0,i1,j0,j1)
1711 : use filter_module,only: filter2
1712 : use edyn_mpi ,only: mp_gatherlons_f3d,mytidi
1713 : use edyn_mpi ,only: mp_scatterlons_f3d
1714 : use edyn_geogrid ,only: nlon
1715 : !
1716 : ! Args:
1717 : integer,intent(in) :: k0,k1,i0,i1,j0,j1
1718 : real(r8),intent(inout) :: f(k0:k1,i0:i1,j0:j1)
1719 : !
1720 : ! Local:
1721 : integer :: i,j,k,nlevs
1722 0 : real(r8) :: fik(nlon,k1-k0+1)
1723 : type(array_ptr_type) :: fkij(1) ! fkij(1)%ptr(k1-k0+1,nlon,j0:j1)
1724 :
1725 0 : nlevs = k1-k0+1
1726 : !
1727 : ! Define lons in fkij from current task subdomain:
1728 : !
1729 0 : allocate(fkij(1)%ptr(nlevs,nlon,j0:j1))
1730 : !$omp parallel do private( i,j,k )
1731 0 : do j=j0,j1
1732 0 : do i=i0,i1
1733 0 : do k=k0,k1
1734 0 : fkij(1)%ptr(k-k0+1,i,j) = f(k,i,j)
1735 : enddo
1736 : enddo
1737 : enddo
1738 : !
1739 : ! Gather longitudes into tasks in first longitude column of task table
1740 : ! (leftmost of each j-row) for global fft. (i.e., tasks with mytidi==0
1741 : ! gather lons from other tasks in that row). This includes all latitudes.
1742 : !
1743 0 : call mp_gatherlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1)
1744 : !
1745 : ! Only leftmost tasks at each j-row of tasks does the global filtering:
1746 : !
1747 0 : if (mytidi==0) then
1748 : !
1749 : ! Define 2d array with all longitudes for filter at each latitude:
1750 : !
1751 0 : do j=j0,j1
1752 0 : do i=1,nlon
1753 0 : do k=k0,k1
1754 0 : fik(i,k-k0+1) = fkij(1)%ptr(k-k0+1,i,j)
1755 : enddo
1756 : enddo
1757 : !
1758 : ! Remove wave numbers > kut(lat):
1759 : !
1760 0 : call filter2(fik,1,nlevs,j)
1761 : !
1762 : ! Return filtered array to fkij:
1763 : !
1764 0 : do i=1,nlon
1765 0 : do k=k0,k1
1766 0 : fkij(1)%ptr(k-k0+1,i,j) = fik(i,k-k0+1)
1767 : enddo
1768 : enddo ! i=1,nlon
1769 : enddo ! j=j0,j1
1770 : endif ! mytidi==0
1771 : !
1772 : ! Now leftmost task at each j-row must redistribute filtered data
1773 : ! back to other tasks in the j-row (mytidi>0,mytidj) (includes latitude):
1774 : !
1775 0 : call mp_scatterlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1)
1776 : !
1777 : ! Return filtered array to inout field at task subdomain:
1778 0 : do j=j0,j1
1779 0 : do i=i0,i1
1780 0 : do k=k0,k1
1781 0 : f(k,i,j) = fkij(1)%ptr(k-k0+1,i,j)
1782 : enddo
1783 : enddo
1784 : enddo
1785 0 : deallocate(fkij(1)%ptr)
1786 0 : end subroutine filter2_op
1787 : !-----------------------------------------------------------------------
1788 : !-----------------------------------------------------------------------
1789 36480 : subroutine ringfilter_op(f,k0,k1,i0,i1,j0,j1)
1790 : use filter_module,only: ringfilter
1791 : use edyn_mpi ,only: mp_gatherlons_f3d,mytidi
1792 : use edyn_mpi ,only: mp_scatterlons_f3d
1793 : use edyn_geogrid ,only: nlon
1794 : !
1795 : ! Args:
1796 : integer,intent(in) :: k0,k1,i0,i1,j0,j1
1797 : real(r8),intent(inout) :: f(k0:k1,i0:i1,j0:j1)
1798 : !
1799 : ! Local:
1800 : integer :: i,j,k,nlevs
1801 72960 : real(r8) :: fik(nlon,k1-k0+1)
1802 : type(array_ptr_type) :: fkij(1) ! fkij(1)%ptr(k1-k0+1,nlon,j0:j1)
1803 :
1804 36480 : nlevs = k1-k0+1
1805 : !
1806 : ! Define lons in fkij from current task subdomain:
1807 : !
1808 182400 : allocate(fkij(1)%ptr(nlevs,nlon,j0:j1))
1809 : !$omp parallel do private( i,j,k )
1810 145920 : do j=j0,j1
1811 1459200 : do i=i0,i1
1812 61833600 : do k=k0,k1
1813 61724160 : fkij(1)%ptr(k-k0+1,i,j) = f(k,i,j)
1814 : enddo
1815 : enddo
1816 : enddo
1817 : !
1818 : ! Gather longitudes into tasks in first longitude column of task table
1819 : ! (leftmost of each j-row) for global fft. (i.e., tasks with mytidi==0
1820 : ! gather lons from other tasks in that row). This includes all latitudes.
1821 : !
1822 36480 : call mp_gatherlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1)
1823 : !
1824 : ! Only leftmost tasks at each j-row of tasks does the global filtering:
1825 : !
1826 36480 : if (mytidi==0) then
1827 : !
1828 : ! Define 2d array with all longitudes for filter at each latitude:
1829 : !
1830 12160 : do j=j0,j1
1831 1322400 : do i=1,nlon
1832 61733280 : do k=k0,k1
1833 61724160 : fik(i,k-k0+1) = fkij(1)%ptr(k-k0+1,i,j)
1834 : enddo
1835 : enddo
1836 : !
1837 : ! Remove wave numbers > kut(lat):
1838 : !
1839 9120 : call ringfilter(1,nlevs,j,fik)
1840 : !
1841 : ! Return filtered array to fkij:
1842 : !
1843 1325440 : do i=1,nlon
1844 61733280 : do k=k0,k1
1845 61724160 : fkij(1)%ptr(k-k0+1,i,j) = fik(i,k-k0+1)
1846 : enddo
1847 : enddo ! i=1,nlon
1848 : enddo ! j=j0,j1
1849 : endif ! mytidi==0
1850 : !
1851 : ! Now leftmost task at each j-row must redistribute filtered data
1852 : ! back to other tasks in the j-row (mytidi>0,mytidj) (includes latitude):
1853 : !
1854 36480 : call mp_scatterlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1)
1855 : !
1856 : ! Return filtered array to inout field at task subdomain:
1857 145920 : do j=j0,j1
1858 1459200 : do i=i0,i1
1859 61833600 : do k=k0,k1
1860 61724160 : f(k,i,j) = fkij(1)%ptr(k-k0+1,i,j)
1861 : enddo
1862 : enddo
1863 : enddo
1864 36480 : deallocate(fkij(1)%ptr)
1865 36480 : end subroutine ringfilter_op
1866 : !-----------------------------------------------------------------------
1867 : end module oplus
|
