module oplus ! ! Horizontally transport the O+ ion, adapted for WACCM-X from TIEGCM. ! Input O+ is received from WACCM physics/chemistry, transported O+ ! (op_out and opnm_out) are passed back to chemistry. ! ! B. Foster (foster@ucar.edu), May, 2015. ! use shr_kind_mod ,only: r8 => shr_kind_r8 use cam_logfile ,only: iulog use pmgrid ,only: plev,plat,plon ! global dimensions use savefield_waccm,only: savefld_waccm, & ! save field to waccm history savefld_waccm_switch use edyn_geogrid ,only: dphi,dlamda,cs,zp,expz,p0 use getapex ,only: bx,by,bz,bmod2 ! (0:nlonp1,jspole-1:jnpole+1) use edyn_params ,only: re use time_manager ,only: get_nstep,get_step_size,is_first_restart_step use cam_abortutils ,only: endrun use edyn_mpi ,only: array_ptr_type use shr_const_mod ,only : shr_const_g ! gravitational constant (m/s^2) implicit none save real(r8) :: pi,rtd ! ! Constants in CGS: ! ! (using shr_const_g for gravitational constant, see use-association above) ! real(r8),parameter :: grav_cm = 870._r8 ! gravity (cm/s2) as in TIEGCM (about 300 km) ! real(r8),parameter :: grav_cm = 980._r8 ! gravity (cm/s2) as in WACCM (surface) real(r8),parameter :: boltz = 1.38E-16_r8 ! boltzman's constant (erg/kelvin) real(r8),parameter :: gask = 8.314e7_r8 ! gas constant (erg/mol) ! ! Collision factor (tuneable) (see also local colfac in iondrag.F90) ! FIX: Collision factor colfac is set locally in iondrag.F90 and here. ! It should be in one location and shared between ionosphere and ! dpie_coupling. ! real(r8),parameter :: colfac = 1.5_r8 ! see also iondrag.F90 ! ! Reciprocal of molecular mass (multiply is cheaper than divide) real(r8),parameter :: rmassinv_o2=1._r8/32._r8, rmassinv_o1=1._r8/16._r8, & rmassinv_n2=1._r8/28._r8 real(r8),parameter :: rmass_op=16._r8 real(r8),parameter :: zbot = 90._r8 ! approx height at bottom of O+ transport (km) real(r8) :: dzp ! delta zp (typically 0.5 from kbot to top) real(r8) :: grav_cm ! gravitational constant (cm/s^2) integer,save :: kbot ! k-index corresponding to ~zbot integer :: nstep ! current model timestep ! ! The shapiro constant .03 is used for spatial smoothing of oplus, ! (shapiro is tuneable, and maybe should be a function of timestep size). ! dtsmooth and dtsmooth_div2 are used in the time smoothing. ! To turn off all smoothing here, set shapiro=0. and dtsmooth = 1. ! real(r8),parameter :: & shapiro = 3.0e-2_r8 ! shapiro constant for spatial smoother real(r8),parameter :: & dtsmooth = 0.95_r8, & ! for time smoother dtsmooth_div2 = 0.5_r8*(1._r8-dtsmooth) ! real(r8),parameter :: & ! shapiro = 0._r8 ! shapiro constant for spatial smoother (0 to turn off) ! real(r8),parameter :: & ! for time smoother ! dtsmooth = 1.0_r8, & ! set to 1.0 to turn off ! dtsmooth_div2 = 0.5_r8*(1._r8-dtsmooth) contains !----------------------------------------------------------------------- subroutine oplus_xport(tnx,tex,tix,unx,vnx,omx,zgx,o2x,o1x,n2x,opx,& opnmx,mbarx,uix,vix,wix,pmid,op_out,opnm_out,& i0,i1,j0,j1,nspltop,ispltop) ! ! All input fields from dpie_coupling are in "TIEGCM" format, i.e., ! longitude (-180->180), vertical (bot2top), and units (CGS). ! use edyn_mpi,only: mp_geo_halos,mp_pole_halos,setpoles use edynamo ,only: find_kbotdyn use filter_module,only: kut2 ! for print only ! ! Transport O+ ion. ! March-May, 2015 B.Foster: Adapted from TIEGCM (oplus.F) for WACCM-X. ! ! Notes: ! - waccmx_opt='ionosphere' must be set in user_nl_cam for te,ti inputs to have values ! ! Args: ! integer,intent(in) :: & i0, & ! grid%ifirstxy i1, & ! grid%ilastxy j0, & ! grid%jfirstxy j1 ! grid%jlastxy integer,intent(in) :: nspltop,ispltop ! ! Input fields without halo points (lon +/-180, vertical bot2top, CGS units): ! real(r8),intent(in) :: tnx (plev,i0:i1,j0:j1) ! neutral temperature (deg K) real(r8),intent(in) :: tex (plev,i0:i1,j0:j1) ! electron temperature (deg K) real(r8),intent(in) :: tix (plev,i0:i1,j0:j1) ! ion temperature (deg K) real(r8),intent(in) :: unx (plev,i0:i1,j0:j1) ! neutral zonal wind (cm/s) real(r8),intent(in) :: vnx (plev,i0:i1,j0:j1) ! neutral meridional wind (cm/s) real(r8),intent(in) :: omx (plev,i0:i1,j0:j1) ! omega (1/s) real(r8),intent(in) :: zgx (plev,i0:i1,j0:j1) ! geopotential height (cm) real(r8),intent(in) :: o2x (plev,i0:i1,j0:j1) ! o2 (mmr) real(r8),intent(in) :: o1x (plev,i0:i1,j0:j1) ! o (mmr) real(r8),intent(in) :: n2x (plev,i0:i1,j0:j1) ! n2 (mmr) real(r8),intent(in) :: opx (plev,i0:i1,j0:j1) ! O+ density (cm^3) real(r8),intent(in) :: opnmx(plev,i0:i1,j0:j1) ! O+ density (cm^3) at time-1 real(r8),intent(in) :: mbarx(plev,i0:i1,j0:j1) ! mean molecular weight ! ! Ion drifts from edynamo (also in tiegcm-format): ! real(r8),intent(in) :: uix(plev,i0:i1,j0:j1) ! zonal ion drift real(r8),intent(in) :: vix(plev,i0:i1,j0:j1) ! meridional ion drift real(r8),intent(in) :: wix(plev,i0:i1,j0:j1) ! vertical ion drift real(r8),intent(in) :: pmid(plev) ! pressure at midpoints (Pa) ! ! Output: ! real(r8),intent(out) :: & op_out (plev,i0:i1,j0:j1), & ! O+ output opnm_out(plev,i0:i1,j0:j1) ! O+ output at time n-1 ! ! Local: ! integer :: i,j,k,lat,k0,k1,jm1,jp1,jm2,jp2,lat0,lat1,kk real(r8),dimension(i0:i1,j0:j1) :: & opflux, & ! upward number flux of O+ (returned by sub oplus_flux) dvb ! divergence of B-field ! ! Local inputs with added halo points in lat,lon: ! real(r8),dimension(plev,i0-2:i1+2,j0-2:j1+2),target :: & tn,te,ti,un,vn,om,zg,op,opnm,o2,o1,n2,mbar real(r8),dimension(plev,i0:i1,j0:j1) :: ui,vi,wi real(r8),dimension(plev,i0-2:i1+2,j0-2:j1+2),target :: & tr ,& ! Reduced temperature (.5*(tn+ti)) tp ,& ! Plasma temperature N(O+)*(te+ti) dj ,& ! diffusion coefficients bvel ,& ! bvel @ j = (B.U)*N(O+) diffj ,& ! (D/(H*DZ)*2.*TP+M*G/R)*N(O+) bdotdh_op ,& ! (b(h)*del(h))*phi bdotdh_opj ,& ! (b(h)*del(h))*phi bdotdh_diff ,& ! (b(h)*del(h))*phi opnm_smooth ! O+ at time-1, smoothed real(r8),dimension(plev,i0:i1,j0:j1) :: & ! for saving to histories diag0,diag1,diag2,diag3,diag4,diag5,diag6,diag7,diag8,diag9,& diag10,diag11,diag12,diag13,diag14,diag15,diag16,diag17,& diag18,diag19,diag20,diag21,diag22,diag23,diag24,diag25,& diag26,diag27 real(r8),dimension(plev,i0:i1,j0-1:j1+1) :: hj ! scale height real(r8) :: gmr,dtime,dtx2,dtx2inv real(r8),dimension(plev,i0:i1) :: & bdzdvb_op, & hdz, & tp1, & tphdz0, & tphdz1, & djint, & divbz, & hdzmbz, & hdzpbz, & bdotu, & tmp1_ki,tmp2_ki ! ! Arguments for tridiagonal solver trsolv (no halos): real(r8),dimension(plev,i0:i1,j0:j1) :: & explicit,explicit_a,explicit_b,p_coeff,q_coeff,r_coeff real(r8),dimension(i0:i1) :: ubca, ubcb ! O+ upper boundary real(r8) :: exmin,exmax real(r8) :: fpole_jpm2(plev,plon,4) real(r8),parameter :: one=1._r8 integer :: kprint logical :: calltrsolv ! ! Pointers for multiple-field calls (e.g., mp_geo_halos) integer,parameter :: nfields=12 real(r8) :: polesign(nfields) type(array_ptr_type),allocatable :: ptrs(:) ! ! Execute: ! nstep = get_nstep() dtime = get_step_size() ! step size in seconds dtime = dtime / float(nspltop) dtx2 = 2._r8*dtime dtx2inv = 1._r8/dtx2 if (nstep==1.and.ispltop==1) then write(iulog,"('oplus: shapiro=',es12.4,' dtsmooth=',es12.4,' dtsmooth_div2=',es12.4)") & shapiro,dtsmooth,dtsmooth_div2 write(iulog,"('oplus: kut2=',/,(20i3))") kut2 write(iulog,"('oplus: shr_const_g=',f8.3)") shr_const_g endif tn=0._r8 ; ti=0._r8 ; te=0._r8 ; un=0._r8 ; vn=0._r8 ; om=0._r8 o2=0._r8 ; o1=0._r8 ; n2=0._r8 ; mbar=0._r8 ; zg=0._r8 ; op=0._r8 opnm=0._r8 ! ! zp,expz are declared in edyn_geogrid.F90, and allocated in sub ! set_geogrid (edyn_init.F90). pmid was passed in here (bot2top) ! from dpie_coupling. ! do k=1,plev zp(k) = -log(pmid(k)*10._r8/p0) expz(k) = exp(-zp(k)) enddo dzp = zp(plev)-zp(plev-1) ! use top 2 levels (typically dzp=0.5) ! write(iulog,"('oplus: plev=',i3,' zp (bot2top) =',/,(6es12.3))") plev,zp ! write(iulog,"('oplus: plev=',i3,' expz (bot2top) =',/,(6es12.3))") plev,expz ! write(iulog,"('oplus: plev=',i3,' dzp =',/,(6es12.3))") plev,dzp ! ! Set subdomain blocks from input (composition is in mmr): ! do k=1,plev do j=j0,j1 do i=i0,i1 tn(k,i,j) = tnx(k,i,j) te(k,i,j) = tex(k,i,j) ti(k,i,j) = tix(k,i,j) un(k,i,j) = unx(k,i,j) vn(k,i,j) = vnx(k,i,j) om(k,i,j) = omx(k,i,j) zg(k,i,j) = zgx(k,i,j) op(k,i,j) = opx(k,i,j) opnm(k,i,j) = opnmx(k,i,j) o2(k,i,j) = o2x(k,i,j) o1(k,i,j) = o1x(k,i,j) n2(k,i,j) = n2x(k,i,j) mbar(k,i,j) = mbarx(k,i,j) ui(k,i,j) = uix(k,i,j) vi(k,i,j) = vix(k,i,j) wi(k,i,j) = wix(k,i,j) enddo enddo enddo ! ! Define halo points on inputs: ! WACCM has global longitude values at the poles (j=1,j=plev) ! (they are constant for most, except the winds.) ! ! Set two halo points in lat,lon: ! real(r8),dimension(plev,i0-2:i1+2,j0-2:j1+2),target :: tn,te,etc. ! allocate(ptrs(nfields)) ! nfields==12 ptrs(1)%ptr => tn ; ptrs(2)%ptr => te ; ptrs(3)%ptr => ti ptrs(4)%ptr => un ; ptrs(5)%ptr => vn ; ptrs(6)%ptr => om ptrs(7)%ptr => o2 ; ptrs(8)%ptr => o1 ; ptrs(9)%ptr => n2 ptrs(10)%ptr => op ; ptrs(11)%ptr => opnm ; ptrs(12)%ptr => mbar polesign = 1._r8 polesign(4:5) = -1._r8 ! un,vn ! ! mp_geo_halos first arg: ! type(array_ptr_type) :: fmsub(nf) ! (lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2) ! call mp_geo_halos(ptrs,1,plev,i0,i1,j0,j1,nfields,polesign) ! ! Set latitude halo points over the poles (this does not change the poles). ! (the 2nd halo over the poles will not actually be used (assuming lat loops ! are lat=2,plat-1), because jp1,jm1 will be the pole itself, and jp2,jm2 ! will be the first halo over the pole) ! ! mp_pole_halos first arg: ! type(array_ptr_type) :: f(nf) ! (plev,i0-2:i1+2,j0-2:j1+2) call mp_pole_halos(ptrs,1,plev,i0,i1,j0,j1,nfields,polesign) deallocate(ptrs) ! ! Use below to exclude the poles (lat=2,plat-1) from latitude scans. ! lat0 = j0 lat1 = j1 if (j0 == 1) lat0 = 2 if (j1 == plat) lat1 = plat-1 ! ! kbot is the k-index at the bottom of O+ transport calculations, ! corresponding to height zbot. kbot is the same for all processors. ! Function find_kbotdyn expects z to be "bot2top" (cm), and will ! search from the bottom at k=1 (function is in the edynamo module) ! Note that kbot is saved, so find_kbotdyn needs to be called only once per run. ! if ((nstep==1.or.is_first_restart_step()).and.ispltop==1) then kbot = find_kbotdyn(zg(:,i0:i1,j0:j1),plev,i0,i1,j0,j1,zbot) write(iulog,"('oplus: zbot=',f8.2,' kbot=',i4,' pmid(kbot)=',es12.4,' zp(kbot)=',es12.4)") & zbot,kbot,pmid(kbot),zp(kbot) endif ! ! Save input fields to WACCM histories. Sub savefld_waccm_switch converts ! fields from tiegcm-format to waccm-format before saving to waccm histories. ! ! call savefld_waccm_switch(tn(:,i0:i1,j0:j1),'OPLUS_TN',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(te(:,i0:i1,j0:j1),'OPLUS_TE',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(ti(:,i0:i1,j0:j1),'OPLUS_TI',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(un(:,i0:i1,j0:j1),'OPLUS_UN',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(vn(:,i0:i1,j0:j1),'OPLUS_VN',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(om(:,i0:i1,j0:j1),'OPLUS_OM',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(zg(:,i0:i1,j0:j1),'OPLUS_Z' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(o2(:,i0:i1,j0:j1),'OPLUS_O2',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(o1(:,i0:i1,j0:j1),'OPLUS_O1',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(n2(:,i0:i1,j0:j1),'OPLUS_N2',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(op(:,i0:i1,j0:j1),'OPLUS_OP',plev,i0,i1,j0,j1) call savefld_waccm_switch(ui(:,i0:i1,j0:j1),'OPLUS_UI',plev,i0,i1,j0,j1) call savefld_waccm_switch(vi(:,i0:i1,j0:j1),'OPLUS_VI',plev,i0,i1,j0,j1) call savefld_waccm_switch(wi(:,i0:i1,j0:j1),'OPLUS_WI',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(mbar(:,i0:i1,j0:j1),'OPLUS_MBAR',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(opnm(:,i0:i1,j0:j1),'OPLUS_OPNM',plev,i0,i1,j0,j1) ! ! Initialize output op_out with input op at 1:kbot-1, to retain values from ! bottom of column up to kbot. This routine will change (transport) these ! outputs only from kbot to the top (plev). ! op_out = 0._r8 opnm_out = 0._r8 op_out (1:kbot-1,i0:i1,j0:j1) = op (1:kbot-1,i0:i1,j0:j1) opnm_out(1:kbot-1,i0:i1,j0:j1) = opnm(1:kbot-1,i0:i1,j0:j1) ! ! Sub oplus_flux returns upward number flux of O+ in opflux ! Output opflux(i,j) is 2d lon x lat subdomain: ! call oplus_flux(opflux,i0,i1,j0,j1) call savefld_waccm(opflux(i0:i1,j0:j1),'OPLUS_FLUX',1,i0,i1,j0,j1) ! ! Divergence of B (mag field) is returned by divb in dvb(i0:i1,j0:j1) ! call divb(dvb,i0,i1,j0,j1) call savefld_waccm(dvb(i0:i1,j0:j1),'OPLUS_DIVB',1,i0,i1,j0,j1) ! ! The solver will be called only if calltrsolv=true. It is sometimes ! set false when skipping parts of the code for debug purposes. ! calltrsolv = .true. tr = 0._r8 tp = 0._r8 dj = 0._r8 hj = 0._r8 bvel= 0._r8 diffj = 0._r8 opnm_smooth = 0._r8 diag0 =0._r8 ! !----------------------- Begin first latitude scan --------------------- do lat=lat0,lat1 jm2 = lat-2 jm1 = lat-1 jp1 = lat+1 jp2 = lat+2 ! write(iulog,"('first lat scan: lat0,1=',2i4,' lat=',i4)") lat0,lat1,lat ! ! as of April, 2015, TIEGCM incorrectly uses te+ti instead of tn+ti ! This has not been fixed in TIEGCM, because fixing it causes a tuning ! problem (ask Hanli and Wenbin). For WACCM, it is correct as below. ! (see also tp) ! do i=i0,i1 ! ! Reduced temperature (tpj in tiegcm): ! 'OPLUS_TR' (has constants at poles) ! do k=kbot,plev tr(k,i,jm1) = 0.5_r8*(tn(k,i,jm1)+ti(k,i,jm1)) tr(k,i,lat) = 0.5_r8*(tn(k,i,lat)+ti(k,i,lat)) tr(k,i,jp1) = 0.5_r8*(tn(k,i,jp1)+ti(k,i,jp1)) enddo enddo ! i=i0,i1 ! ! rrk returns ambipolar diffusion coefficients in d(jm1),dj(lat),djp1(jp1): ! 'OPLUS_DJ' (has constants at poles) ! call rrk( & tn(kbot:plev,i0:i1,jm1),mbar(kbot:plev,i0:i1,jm1), & o2(kbot:plev,i0:i1,jm1),o1 (kbot:plev,i0:i1,jm1), & n2(kbot:plev,i0:i1,jm1),tr (kbot:plev,i0:i1,jm1), & dj(kbot:plev,i0:i1,jm1),i0,i1,kbot,plev,lat) call rrk( & tn(kbot:plev,i0:i1,lat),mbar(kbot:plev,i0:i1,lat), & o2(kbot:plev,i0:i1,lat),o1 (kbot:plev,i0:i1,lat), & n2(kbot:plev,i0:i1,lat),tr (kbot:plev,i0:i1,lat), & dj(kbot:plev,i0:i1,lat),i0,i1,kbot,plev,lat) call rrk( & tn(kbot:plev,i0:i1,jp1),mbar(kbot:plev,i0:i1,jp1), & o2(kbot:plev,i0:i1,jp1),o1 (kbot:plev,i0:i1,jp1), & n2(kbot:plev,i0:i1,jp1),tr (kbot:plev,i0:i1,jp1), & dj(kbot:plev,i0:i1,jp1),i0,i1,kbot,plev,lat) ! ! Plasma temperature: ! 'OPLUS_TP0' (tp will get poles from jm1 and jp1) ! do i=i0,i1 do k=kbot,plev tp(k,i,jm1) = te(k,i,jm1)+ti(k,i,jm1) tp(k,i,lat) = te(k,i,lat)+ti(k,i,lat) tp(k,i,jp1) = te(k,i,jp1)+ti(k,i,jp1) enddo enddo diag0(kbot:plev,i0:i1,lat) = tp(kbot:plev,i0:i1,lat) ! ! Add poles to diag0: if (j0==1.and.lat==2) diag0(kbot:plev,i0:i1,j0) = tp(kbot:plev,i0:i1,jm1) if (j1==plat.and.lat==plat-1) diag0(kbot:plev,i0:i1,j1) = tp(kbot:plev,i0:i1,jp1) ! ! Neutral scale height: ! 'OPLUS_HJ' (has constants at poles) ! grav_cm = shr_const_g * 100._r8 ! m/s^2 -> cm/s^2 do i=i0,i1 do k=kbot,plev hj(k,i,jm1) = gask * tn(k,i,jm1) / (mbar(k,i,jm1) * grav_cm) hj(k,i,lat) = gask * tn(k,i,lat) / (mbar(k,i,lat) * grav_cm) hj(k,i,jp1) = gask * tn(k,i,jp1) / (mbar(k,i,jp1) * grav_cm) enddo enddo ! ! bvel @ jm1 = (B.U)*N(O+) (J-1) ! bvel @ j = (B.U)*N(O+) (J) ! bvel @ jp1 = (B.U)*N(O+) (J+1) ! 'OPLUS_BVEL' (has constants at poles) ! ! Note bx,by,bz were set globally for all tasks by sub magfield ! (getapex.F90) ! do i=i0,i1 do k=kbot,plev bvel(k,i,jm1) = & (bx(i,jm1)*un(k,i,jm1)+by(i,jm1)*vn(k,i,jm1)+ & hj(k,i,jm1)*bz(i,jm1)*om(k,i,jm1))*op(k,i,jm1) bvel(k,i,lat) = & (bx(i,lat)*un(k,i,lat)+by(i,lat)*vn(k,i,lat)+ & hj(k,i,lat)*bz(i,lat)*om(k,i,lat))*op(k,i,lat) bvel(k,i,jp1) = & (bx(i,jp1)*un(k,i,jp1)+by(i,jp1)*vn(k,i,jp1)+ & hj(k,i,jp1)*bz(i,jp1)*om(k,i,jp1))*op(k,i,jp1) enddo ! k=kbot,plev enddo ! i=lon0,lon1 ! ! Ambipolar diffusion is returned in diffj: ! 'OPLUS_DIFFJ' (will have constants at poles after this lat scan) ! call diffus(tp(kbot:plev,i0:i1,jm1),op(kbot:plev,i0:i1,jm1),hj(kbot:plev,:,jm1), & diffj(kbot:plev,i0:i1,jm1),i0,i1,kbot,plev,lat) call diffus(tp(kbot:plev,i0:i1,lat),op(kbot:plev,i0:i1,lat),hj(kbot:plev,:,lat), & diffj(kbot:plev,i0:i1,lat),i0,i1,kbot,plev,lat) call diffus(tp(kbot:plev,i0:i1,jp1),op(kbot:plev,i0:i1,jp1),hj(kbot:plev,:,jp1), & diffj(kbot:plev,i0:i1,jp1),i0,i1,kbot,plev,lat) ! ! 'OPLUS_TP1' (constants at the poles) ! do i=i0,i1 do k=kbot,plev tp(k,i,jm2) = op(k,i,jm2)*(te(k,i,jm2)+ti(k,i,jm2)) tp(k,i,jm1) = tp(k,i,jm1)*op(k,i,jm1) tp(k,i,lat) = tp(k,i,lat)*op(k,i,lat) tp(k,i,jp1) = tp(k,i,jp1)*op(k,i,jp1) tp(k,i,jp2) = op(k,i,jp2)*(te(k,i,jp2)+ti(k,i,jp2)) enddo enddo ! ! Latidinal shapiro smoother: opnm is O+ at time n-1. ! opnm_smooth will be used in explicit terms below. ! Smooth in latitude: ! 'OPNM_SMOOTH' (zero at poles) ! do i=i0,i1 do k=kbot,plev opnm_smooth(k,i,lat) = opnm(k,i,lat)-shapiro* & (opnm(k,i,jp2)+opnm(k,i,jm2)-4._r8* & (opnm(k,i,jp1)+opnm(k,i,jm1))+6._r8* & opnm(k,i,lat)) enddo ! k=kbot,plev enddo ! i=i0,i1 enddo ! end first latitude scan (lat=lat0,lat1) ! !------------------------- End first latitude scan --------------------- ! ! Set pole values for opnm_smooth. Do this before savefld calls, so plots will ! include the poles. All other fields in 1st lat scan got values at the poles ! via jm1,jp1 above. ! call setpoles(opnm_smooth(kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! ! Save to history file (exclude halo points) ! ! call savefld_waccm_switch(tr (:,i0:i1,j0:j1),'OPLUS_TR' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(dj (:,i0:i1,j0:j1),'OPLUS_DJ' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(hj (:,i0:i1,j0:j1),'OPLUS_HJ' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(bvel (:,i0:i1,j0:j1),'OPLUS_BVEL' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diffj(:,i0:i1,j0:j1),'OPLUS_DIFFJ',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag0(:,i0:i1,j0:j1),'OPLUS_TP0' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(tp (:,i0:i1,j0:j1),'OPLUS_TP1' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(opnm_smooth(:,i0:i1,j0:j1),'OPNM_SMOOTH',plev,i0,i1,j0,j1) ! ! Set halo points where needed. ! allocate(ptrs(5)) ptrs(1)%ptr => dj ; ptrs(2)%ptr => bvel ; ptrs(3)%ptr => diffj ptrs(4)%ptr => tp ; ptrs(5)%ptr => opnm_smooth polesign(1:5) = 1._r8 call mp_geo_halos (ptrs,1,plev,i0,i1,j0,j1,5,polesign(1:5)) call mp_pole_halos(ptrs,1,plev,i0,i1,j0,j1,5,polesign(1:5)) deallocate(ptrs) ! ! For debug: return after first lat scan: ! The above fields look good if we return here. ! write(iulog,"('oplus returning after first lat scan')") ! return ! !----------------------- Begin second latitude scan -------------------- bdotdh_op = 0._r8 bdotdh_opj = 0._r8 do lat=lat0,lat1 jm2 = lat-2 jm1 = lat-1 jp1 = lat+1 jp2 = lat+2 ! write(iulog,"('second lat scan: lat0,1=',2i4,' lat=',i4)") lat0,lat1,lat ! ! bdotdh_op = (B(H).DEL(H))*(D/(H*DZ)*TP+M*G/R)*N(O+) ! then bdotdh_op = d*bz*bdotdh_op ! real(r8),dimension(plev,i0-2:i1+2,j0-2:j1+2) :: diffj ! real(r8),dimension(plev,i0-2:i1+2,j0-2:j1+2) :: bdotdh_op ! 'BDOTDH_OP' (zero at the poles) ! call bdotdh( & diffj(kbot:plev,i0:i1,jm1), & diffj(kbot:plev,:,lat ), & ! includes longitude halos diffj(kbot:plev,i0:i1,jp1), & bdotdh_op(kbot:plev,i0:i1,lat),i0,i1,kbot,plev,lat) ! do i=i0,i1 do k=kbot,plev bdotdh_op(k,i,lat) = dj(k,i,lat)*bz(i,lat)*bdotdh_op(k,i,lat) ! BDOTDH_OP enddo ! k=kbot,plev enddo ! i=i0,i1 ! ! bdotdh_opjm1 = (B(H).DEL(H))*2.*TP*N(O+) (J-1) ! bdotdh_opj = (B(H).DEL(H))*2.*TP*N(O+) (J) ! bdotdh_opjp1 = (B(H).DEL(H))*2.*TP*N(O+) (J+1) ! 'BDOTDH_OPJ' (has reasonable non-constant values at poles) ! call bdotdh( & tp(kbot:plev,i0:i1,jm2), & tp(kbot:plev,:,jm1), & tp(kbot:plev,i0:i1,lat), & bdotdh_opj(kbot:plev,i0:i1,jm1),i0,i1,kbot,plev,jm1) call bdotdh( & tp(kbot:plev,i0:i1,jm1), & tp(kbot:plev,:,lat), & tp(kbot:plev,i0:i1,jp1), & bdotdh_opj(kbot:plev,i0:i1,lat),i0,i1,kbot,plev,lat) call bdotdh( & tp(kbot:plev,i0:i1,lat), & tp(kbot:plev,:,jp1), & tp(kbot:plev,i0:i1,jp2), & bdotdh_opj(kbot:plev,i0:i1,jp1),i0,i1,kbot,plev,jp1) ! do i=i0,i1 do k=kbot,plev bdotdh_opj(k,i,jm1) = bdotdh_opj(k,i,jm1)*dj(k,i,jm1) bdotdh_opj(k,i,lat) = bdotdh_opj(k,i,lat)*dj(k,i,lat) bdotdh_opj(k,i,jp1) = bdotdh_opj(k,i,jp1)*dj(k,i,jp1) enddo ! k=kbot,plev enddo ! i=i0,i1 enddo ! lat=j0,j1 (end second lat scan) ! !------------------------ End second latitude scan --------------------- ! ! bdotdh_opj already has non-constant polar values, but bdotdh_op poles are zero. ! Sub setpoles will set poles to the zonal average of the latitude below each pole. ! ! This may not be necessary, but do it for plotting: call setpoles(bdotdh_op(kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! kprint = kbot+5 ! call printpoles(bdotdh_op (kbot:plev,i0:i1,j0:j1),kprint,kbot,plev,i0,i1,j0,j1,'BDOTDH_OP' ) ! call printpoles(bdotdh_opj(kbot:plev,i0:i1,j0:j1),kprint,kbot,plev,i0,i1,j0,j1,'BDOTDH_OPJ') ! call savefld_waccm_switch(bdotdh_op (:,i0:i1,j0:j1),'BDOTDH_OP' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(bdotdh_opj(:,i0:i1,j0:j1),'BDOTDH_OPJ',plev,i0,i1,j0,j1) ! ! For debug: return after second lat scan: ! write(iulog,"('oplus returning after second lat scan')") ! return ! ! Note mp_geo_halos will overwrite jm1,jp1 that was set above. ! bdotdh_opj needs longitude halos for the bdotdh call below. ! ! real(r8),dimension(plev,i0-2:i1+2,j0-2:j1+2),target :: bdotdh_op,opj ! allocate(ptrs(1)) ptrs(1)%ptr => bdotdh_opj call mp_geo_halos (ptrs,1,plev,i0,i1,j0,j1,1,(/1._r8/)) call mp_pole_halos(ptrs,1,plev,i0,i1,j0,j1,1,(/1._r8/)) deallocate(ptrs) ! !----------------------- Begin third latitude scan --------------------- ! bdotdh_diff = 0._r8 bdzdvb_op = 0._r8 explicit = 0._r8 ; explicit_a=0._r8 ; explicit_b=0._r8 hdz = 0._r8 tphdz0 = 0._r8 tphdz1 = 0._r8 djint = 0._r8 divbz = 0._r8 hdzmbz = 0._r8 hdzpbz = 0._r8 p_coeff = 0._r8 q_coeff = 0._r8 r_coeff = 0._r8 bdotu = 0._r8 diag1 = 0._r8 ; diag2 = 0._r8 ; diag3 = 0._r8 ; diag4 = 0._r8 ; diag5 = 0._r8 diag6 = 0._r8 ; diag7 = 0._r8 ; diag8 = 0._r8 ; diag9 = 0._r8 ; diag10= 0._r8 diag11 = 0._r8 ; diag12= 0._r8 ; diag13= 0._r8 ; diag14= 0._r8 ; diag15= 0._r8 diag16 = 0._r8 ; diag17= 0._r8 ; diag18= 0._r8 ; diag19= 0._r8 ; diag20= 0._r8 diag21 = 0._r8 ; diag22= 0._r8 ; diag23= 0._r8 ; diag24= 0._r8 ; diag25= 0._r8 diag26 = 0._r8 ; diag27= 0._r8 ! ! Globally, this loop is lat=2,plat-1 (i.e., skipping the poles) ! do lat=lat0,lat1 jm2 = lat-2 jm1 = lat-1 ! this will be south pole for southern pes (j==1) jp1 = lat+1 ! this will be north pole for northern pes (j==plat) jp2 = lat+2 ! write(iulog,"('third lat scan: lat0,1=',2i4,' lat=',i4)") lat0,lat1,lat ! ! bdotdh_opj = (B(H).DEL(H))*D*(B(H).DEL(H))*2.*TP*N(O+) (J) ! 'BDOTDH_DIFF' (zero at the poles) ! call bdotdh( & bdotdh_opj(kbot:plev,i0:i1,jm1), & bdotdh_opj(kbot:plev,:,lat), & ! includes longitude halos bdotdh_opj(kbot:plev,i0:i1,jp1), & bdotdh_diff(kbot:plev,i0:i1,lat),i0,i1,kbot,plev,lat) ! BDOTDH_DIFF ! ! bdzdvb_op = (BZ*D/(H*DZ)+DIV(*B))*S2 ! bdzdvb returns bdzdvb_op(k,i). ! 'BDZDVB_OP' (zero at the poles) ! ! real(r8),dimension(i0:i1,j0:j1) :: dvb ! real(r8),dimension(plev,i0:i1,j0-1:j1+1) :: hj ! scale height ! real(r8),dimension(plev,i0:i1) :: bdzdvb_op ! subroutine bdzdvb(phi,dvb,h,ans,lev0,lev1,lon0,lon1,lat) ! real(r8),intent(in) :: dvb(lon0:lon1) ! real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: phi,h ! real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: ans ! call bdzdvb(bdotdh_opj(kbot:plev,i0:i1,lat),dvb(:,lat),hj(kbot:plev,i0:i1,lat), & bdzdvb_op(kbot:plev,i0:i1),kbot,plev,i0,i1,lat) diag1(:,i0:i1,lat) = bdzdvb_op(:,i0:i1) ! BDZDVB_OP ! ! Collect explicit terms: ! 'EXPLICIT0' (this will have poles set after third lat scan, before ! plotting. The poles will be constant in longitude, and ! may differ structurally from adjacent latitudes. ! do i=i0,i1 do k=kbot,plev explicit(k,i,lat) = -one*(bdzdvb_op(k,i)+bdotdh_diff(k,i,lat)+ & bdotdh_op(k,i,lat)) enddo ! k=kbot,plev enddo ! i=i0,i1 diag2(:,i0:i1,lat) = explicit(:,i0:i1,lat) ! EXPLICIT0 ! ! For debug: skip rest of this lat loop: ! write(iulog,"('cycling third lat loop after EXPLICIT0')") ! calltrsolv = .false. ! goto 300 ! ! Ion drifts are interpolated to midpoints (is this necessary in WACCM?). ! ! Need lon,lat halos for op, bvel, and bmod2 ! op,bvel halos were set above, bmod2 was set in magfield (getapex.F90) ! (ui,vi,wi halos are not used here.) ! ! bmod2 halos are set in sub magfield (getapex.F90), including plat-1,plat,plat+1, ! and 1 halo point in longitude. Note bmod2 is global in lon and lat for all pe's. ! use getapex,only: bmod2 ! (0:nlonp1,jspole-1:jnpole+1) ! ! When looping lat=2,plat-1, this explicit has zero pole values, ! but there are still problems at processor longitude boundaries, ! especially near the south pole: ! 'EXPLICIT1' (zero at the poles) ! do i=i0,i1 do k=kbot,plev-1 ! ! Original TIEGCM statement: ! explicit(k,i) = explicit(k,i)+1._r8/(2._r8*re)* & ! (1._r8/(cs(lat)*dlamda)*(bx(i,lat)* & ! (bvel(k,i+1,lat)-bvel(k,i-1,lat))+ & ! 0.5_r8*(ui(k,i,lat)+ui(k+1,i,lat))*bmod2(i,lat)**2* & ! (op(k,i+1,lat)/bmod2(i+1,lat)**2- & ! op(k,i-1,lat)/bmod2(i-1,lat)**2))+ & ! ! 1._r8/dphi*(by(i,lat)*(bvel(k,i,jp1)-bvel(k,i,jm1))+ & ! 0.5_r8*(vi(k,i,lat)+vi(k+1,i,lat))*bmod2(i,lat)**2* & ! (op(k,i,jp1)/bmod2(i,jp1)**2- & ! op(k,i,jm1)/bmod2(i,jm1)**2))) ! ! Break it into two pieces and put together for debug: ! ! 'EXPLICITa' explicit_a(k,i,lat) = (bx(i,lat)* & (bvel(k,i+1,lat)-bvel(k,i-1,lat))+ & 0.5_r8*(ui(k,i,lat)+ui(k+1,i,lat))*bmod2(i,lat)**2* & (op(k,i+1,lat)/bmod2(i+1,lat)**2- & op(k,i-1,lat)/bmod2(i-1,lat)**2)) ! ! 'EXPLICITb' ! explicit_b(k,i,lat) = & (by(i,lat)*(bvel(k,i,jp1)-bvel(k,i,jm1))+ & 0.5_r8*(vi(k,i,lat)+vi(k+1,i,lat))*bmod2(i,lat)**2* & (op(k,i,jp1)/bmod2(i,jp1)**2- & op(k,i,jm1)/bmod2(i,jm1)**2)) ! ! 'EXPLICIT1' ! explicit will receive polar values after this latitude scan. ! explicit(k,i,lat) = explicit(k,i,lat)+1._r8/(2._r8*re)* & (1._r8/(cs(lat)*dlamda)*explicit_a(k,i,lat)+ & 1._r8/dphi*explicit_b(k,i,lat)) ! ! explicit is bad at i=1,72,73,144 near south pole (npole appears to be ok) ! This does not appear to adversely affect the final O+ output, and TIEGCM ! has the same high magnitudes, so am ignoring this for now. The high magnitudes ! are near the south pole, at processor longitude boundaries (implicating an error ! with longitude halo points). ! ! if (explicit(k,i,lat) < -300._r8 .or. explicit(k,i,lat) > 300._r8) then ! write(iulog,"('>>> bad explicit: nstep=',i3,' k,i,lat=',3i4,' explicit=',es12.4)") & ! nstep,k,i,lat,explicit(k,i) ! write(iulog,"(' cs(lat) =',3es12.4)") cs(lat) ! write(iulog,"(' op(k,i-1:i+1,lat) =',3es12.4)") op(k,i-1:i+1,lat) ! write(iulog,"(' op(k,i,jm1:jp1) =',3es12.4)") op(k,i,jm1:jp1) ! write(iulog,"(' bvel(k,i-1:i+1,lat)=',3es12.4)") bvel(k,i-1:i+1,lat) ! write(iulog,"(' bvel(k,i,jm1:jp1) =',3es12.4)") bvel(k,i,jm1:jp1) ! write(iulog,"(' bmod2(i-1:i+1,lat) =',3es12.4)") bmod2(i-1:i+1,lat) ! write(iulog,"(' bmod2(i,jm1:jp1) =',3es12.4)") bmod2(i,jm1:jp1) ! write(iulog,"(' ui(k:k+1,i,lat) =',2es12.4)") ui(k:k+1,i,lat) ! write(iulog,"(' vi(k:k+1,i,lat) =',2es12.4)") vi(k:k+1,i,lat) ! write(iulog,"(' bx,by(i,lat) =',2es12.4)") bx(i,lat),by(i,lat) ! endif enddo ! k=kbot,plev-1 enddo ! i=i0,i1 do k=kbot,plev diag25(k,i0:i1,lat) = bmod2(i0:i1,lat) ! BMOD2 (redundant in vertical) enddo diag26(:,i0:i1,lat) = explicit_a(:,i0:i1,lat) ! EXPLICITa diag27(:,i0:i1,lat) = explicit_b(:,i0:i1,lat) ! EXPLICITb diag3 (:,i0:i1,lat) = explicit (:,i0:i1,lat) ! EXPLICIT1 ! ! For debug: skip rest of this lat loop: ! write(iulog,"('cycling third lat loop after EXPLICIT1')") ! goto 300 do i=i0,i1 dvb(i,lat) = dvb(i,lat)/bz(i,lat) enddo ! i=i0,i1 do i=i0,i1 do k=kbot,plev hdz(k,i) = 1._r8/(hj(k,i,lat)*dzp) tp1(k,i) = 0.5_r8*(ti(k,i,lat)+te(k,i,lat)) enddo ! k=kbot,plev enddo ! i=i0,i1 ! ! gmr = G*M(O+)/(2.*R) ! gmr = grav_cm*rmass_op/(2._r8*gask) do i=i0,i1 do k=kbot,plev-1 tphdz1(k+1,i) = 2._r8*tp1(k+1,i)*(0.5_r8*(hdz(k,i)+hdz(k+1,i)))+gmr tphdz0(k+1,i) = 2._r8*tp1(k ,i)*(0.5_r8*(hdz(k,i)+hdz(k+1,i)))-gmr enddo ! k=kbot,plev-1 enddo ! i=lon0,lon1 ! ! Upper and lower boundaries: ! Both TPHDZ0 and TPHDZ1 are zero at the poles. ! ! 5/9/15: Appears to be a problem in TPHDZ0,1 near kbot, maybe is zero? ! do i=i0,i1 tphdz1(kbot,i) = 2._r8*tp1(kbot,i)* & (1.5*hdz(kbot,i)-0.5_r8*hdz(kbot+1,i))+gmr tphdz1(plev,i) = 2._r8*(2._r8*tp1(plev-1,i)-tp1(plev-2,i))* & (1.5_r8*hdz(plev-1,i)-0.5*hdz(plev-2,i))+gmr tphdz0(kbot,i) = 2._r8*(2._r8*tp1(kbot,i)-tp1(kbot+1,i))* & (1.5_r8*hdz(kbot,i)-0.5_r8*hdz(kbot+1,i))-gmr tphdz0(plev,i) = 2._r8*tp1(plev-1,i)* & (1.5_r8*hdz(plev-1,i)-0.5_r8*hdz(plev-2,i))-gmr enddo ! i=i0,i1 diag4(:,i0:i1,lat) = tphdz0(:,i0:i1) ! TPHDZ0 diag5(:,i0:i1,lat) = tphdz1(:,i0:i1) ! TPHDZ1 ! ! djint = dj diffusion at interfaces: ! 'DJINT' (zero at the poles - messes up the plots - may give ! diag6 polar values after the lat scan, before plotting) ! do i=i0,i1 do k=kbot,plev-1 djint(k+1,i) = 0.5_r8*(dj(k,i,lat)+dj(k+1,i,lat)) enddo djint(kbot,i) = (1.5_r8*dj(kbot ,i,lat)-0.5_r8*dj(kbot+1,i,lat)) djint(plev,i) = (1.5_r8*dj(plev-1,i,lat)-0.5_r8*dj(plev-2,i,lat)) enddo ! i=i0,i1 diag6(:,i0:i1,lat) = djint(:,i0:i1) ! DJINT ! ! divbz = (DIV(B)+(DH*D*BZ)/(D*BZ) ! 'DIVBZ' Field appears as a line following mins along magnetic equator (zero at poles) ! Field may be zero at kbot? Lat slices look strange between +/- 14 deg lat. ! do i=i0,i1 do k=kbot,plev divbz(k,i) = & dvb(i,lat)+1._r8/(re*dj(k,i,lat)*bz(i,lat)**2)*(bx(i,lat)/ & cs(lat)*(dj(k,i+1,lat)*bz(i+1,lat)-dj(k,i-1,lat)* & bz(i-1,lat))/(2._r8*dlamda)+by(i,lat)*(dj(k,i,jp1)* & bz(i,jp1)-dj(k,i,jm1)*bz(i,jm1))/(2._r8*dphi)) enddo ! k=kbot,plev enddo ! i=i0,i1 diag7(:,i0:i1,lat) = divbz(:,i0:i1) ! DIVBZ ! ! hdzmbz = (1./(H*DZ)-(DIV(B)+DH*D*BZ/(D*BZ))/(2*BZ))*BZ**2 ! hdzpbz = (1./(H*DZ)+(DIV(B)+DH*D*BZ/(D*BZ))/(2*BZ))*BZ**2 ! 'HDZMBZ' and 'HDZPBZ' are zero at the poles. ! do i=i0,i1 do k=kbot,plev hdzmbz(k,i) = (hdz(k,i)-0.5_r8*divbz(k,i))*bz(i,lat)**2 hdzpbz(k,i) = (hdz(k,i)+0.5_r8*divbz(k,i))*bz(i,lat)**2 enddo ! k=kbot,plev enddo ! i=i0,i1 diag8(:,i0:i1,lat) = hdzmbz(:,i0:i1) ! HDZMBZ diag9(:,i0:i1,lat) = hdzpbz(:,i0:i1) ! HDZPBZ ! ! Sum O+ at time n-1 to explicit terms: N(O+)/(2*DT) (N-1) ! 'EXPLICIT2' (zero at the poles) ! do i=i0,i1 do k=kbot,plev explicit(k,i,lat) = explicit(k,i,lat)-(opnm_smooth(k,i,lat)-shapiro* & (opnm_smooth(k,i+2,lat)+opnm_smooth(k,i-2,lat)-4._r8* & (opnm_smooth(k,i+1,lat)+opnm_smooth(k,i-1,lat))+6._r8* & opnm_smooth(k,i,lat)))*dtx2inv enddo ! k=kbot,plev enddo ! i=i0,i1 diag10(:,i0:i1,lat) = explicit(:,i0:i1,lat) ! EXPLICIT2 ! ! For debug: skip rest of this lat loop: ! write(iulog,"('cycling third lat loop after EXPLICIT2')") ! goto 300 ! ! Begin coefficients p_coeff, q_coeff, r_coeff ! do i=i0,i1 do k=kbot,plev-1 p_coeff(k,i,lat) = hdzmbz(k,i)*djint(k ,i)*tphdz0(k ,i) q_coeff(k,i,lat) = -(hdzpbz(k,i)*djint(k+1,i)*tphdz0(k+1,i)+ & hdzmbz(k,i)*djint(k ,i)*tphdz1(k ,i)) r_coeff(k,i,lat) = hdzpbz(k,i)*djint(k+1,i)*tphdz1(k+1,i) enddo ! k=kbot,plev-1 enddo ! i=i0,i1 diag11(:,i0:i1,lat) = p_coeff(:,i0:i1,lat) ! P_COEFF0 (zero at poles) diag12(:,i0:i1,lat) = q_coeff(:,i0:i1,lat) ! Q_COEFF0 (zero at ubc) diag13(:,i0:i1,lat) = r_coeff(:,i0:i1,lat) ! R_COEFF0 (zero at ubc) ! ! For debug: ! write(iulog,"('cycling third lat loop after P_COEFF0, et.al.')") ! goto 300 ! ! bdotu = B.U ! Introducing neutral winds. ! Am not using 0.5*(om(k)+om(k+1)) here because waccm omega is on midpoints (?) ! (tiegcm has 0.5*(w(k,i,j0)+w(k+1,i,j0)) ) ! do i=i0,i1 do k=kbot,plev bdotu(k,i) = bx(i,lat)*un(k,i,lat)+by(i,lat)*vn(k,i,lat)+ & hj(k,i,lat)*bz(i,lat)*om(k,i,lat) enddo ! k=kbot,plev enddo ! i=i0,i1 diag14(:,i0:i1,lat) = bdotu(:,i0:i1) ! BDOTU ! ! Continue coefficients with vertical ion drift: ! wi is converted from interfaces to midpoints (first use of wi). ! The p,q,r coeffs are still zero at top boundary k=plev, and at poles. ! do i=i0,i1 do k=kbot,plev-2 p_coeff(k+1,i,lat) = p_coeff(k+1,i,lat)+(bz(i,lat)*bdotu(k,i)+ & 0.5_r8*(wi(k+1,i,lat)+wi(k+2,i,lat)))*0.5_r8*hdz(k+1,i) q_coeff(k,i,lat) = q_coeff(k,i,lat)-0.5_r8*(wi(k,i,lat)+wi(k+1,i,lat))*6._r8/re r_coeff(k,i,lat) = r_coeff(k,i,lat)-(bz(i,lat)*bdotu(k+1,i)+ & 0.5_r8*(wi(k,i,lat)+wi(k+1,i,lat)))*0.5_r8*hdz(k,i) enddo ! k=kbot,plev-1 enddo ! i=i0,i1 diag22(:,i0:i1,lat) = p_coeff(:,i0:i1,lat) ! P_COEFF0a diag23(:,i0:i1,lat) = q_coeff(:,i0:i1,lat) ! Q_COEFF0a diag24(:,i0:i1,lat) = r_coeff(:,i0:i1,lat) ! R_COEFF0a ! ! For debug: skip rest of this lat loop: ! p,q,r coeffs look good (zero upper boundaries) if we skip out here: ! write(iulog,"('cycling third lat loop after P_COEFF0a, et.al.')") ! goto 300 ! ! Upper (plev) and lower (kbot) boundaries of p,q,r_coeff: ! (convert wi to midpoints) ! ! tiegcm considers nlev-1 to be the top level. Do it tiegcm-style here, ! and then extrapolate to plev. ! do i=i0,i1 p_coeff(kbot,i,lat) = p_coeff(kbot,i,lat)+(bz(i,lat)* & ! reset p_coeff lbc (2._r8*bdotu(kbot,i)-bdotu(kbot+1,i))+ & 0.5_r8*(wi(kbot,i,lat)+wi(kbot+1,i,lat)))*0.5*hdz(kbot,i) q_coeff(plev-1,i,lat) = q_coeff(plev-1,i,lat)- & 0.5_r8*(wi(plev,i,lat)+wi(plev-1,i,lat))*6._r8/re r_coeff(plev-1,i,lat) = r_coeff(plev-1,i,lat)-(bz(i,lat)* & (2._r8*bdotu(plev-1,i)-bdotu(plev-2,i))+ & 0.5_r8*(wi(plev,i,lat)+wi(plev-1,i,lat)))*0.5_r8*hdz(plev-1,i) enddo ! i=i0,i1 ! ! Extrapolate to top level (tiegcm does not do this): ! p_coeff(plev,i0:i1,lat) = 1.5_r8*p_coeff(plev-1,i0:i1,lat)- & 0.5_r8*p_coeff(plev-2,i0:i1,lat) q_coeff(plev,i0:i1,lat) = 1.5_r8*q_coeff(plev-1,i0:i1,lat)- & 0.5_r8*q_coeff(plev-2,i0:i1,lat) r_coeff(plev,i0:i1,lat) = 1.5_r8*r_coeff(plev-1,i0:i1,lat)- & 0.5_r8*r_coeff(plev-2,i0:i1,lat) ! ! All P,Q,R are zero at the poles. Polar values will be set after third lat scan. diag15(:,i0:i1,lat) = p_coeff(:,i0:i1,lat) ! P_COEFF1 (zero at ubc and poles) diag17(:,i0:i1,lat) = r_coeff(:,i0:i1,lat) ! R_COEFF1 (ok at ubc, zero at poles) ! ! For debug: skip rest of this lat loop: ! write(iulog,"('oplus cycling third lat scan after P_COEFF1, et.al.')") ! goto 300 ! ! Additions to Q coefficients (includes q_coeff lbc,ubc): do i=i0,i1 do k=kbot,plev q_coeff(k,i,lat) = q_coeff(k,i,lat)-bdotu(k,i)*dvb(i,lat)*bz(i,lat)-dtx2inv enddo ! k=kbot,plev-1 enddo ! i=i0,i1 ! ! Plot Q_COEFF1 after ubc has been set. diag16(:,i0:i1,lat) = q_coeff(:,i0:i1,lat) ! Q_COEFF1 (ok at ubc, zero at poles) ! For debug: skip rest of this lat loop: ! write(iulog,"('oplus cycling third lat scan after P_COEFF1, et.al.')") ! goto 300 ! ! Upper boundary condition for O+: do i=i0,i1 ubca(i) = 0._r8 ubcb(i) = -bz(i,lat)**2*djint(plev,i)*tphdz0(plev,i)-ubca(i) ubca(i) = -bz(i,lat)**2*djint(plev,i)*tphdz1(plev,i)+ubca(i) ! ! Q = Q+B/A*R q_coeff(plev,i,lat) = q_coeff(plev,i,lat)+ubcb(i)/ubca(i)* & r_coeff(plev,i,lat) ! ! F = F -R/A*PHI explicit(plev,i,lat) = explicit(plev,i,lat)-opflux(i,lat)* & ! explicit ubc r_coeff(plev,i,lat)/ubca(i) r_coeff(plev,i,lat) = 0._r8 ! r_coeff ubc is reset to zero enddo ! i=i0,i1 ! ! Ubc of EXPLICIT3 has a stripe along the mag equator, unlike the level below. ! diag18(:,i0:i1,lat) = explicit(:,i0:i1,lat) ! EXPLICIT3 (ubc ok, zero at poles) diag19(:,i0:i1,lat) = p_coeff(:,i0:i1,lat) ! P_COEFF2 (zero at ubc, zero at poles) diag20(:,i0:i1,lat) = q_coeff(:,i0:i1,lat) ! Q_COEFF2 (ubc ok, zero at poles) diag21(:,i0:i1,lat) = r_coeff(:,i0:i1,lat) ! R_COEFF2 (zero at ubc, zero at poles) ! ! For debug: skip rest of this lat loop: ! write(iulog,"('cycling third lat loop after P_COEFF2, et.al.')") ! goto 300 ! ! At this point, TIEGCM calculates "sources and sinks" xiop2p and xiop2d. ! Also calculates op_loss, which is subtracted from q_coeff. ! Then TIEGCM "Add source term to RHS (explicit terms)", and calculates ! lower boundary condition N(O+) = Q/L (q_coeff, explicit, p_coeff), and ! finally calls trsolv. ! 300 continue enddo ! end third latitude scan (lat=lat0,lat1) ! ! For debug: return before fourth latitude scan: ! write(iulog,"('oplus returning after third lat scan')") ! op_out (kbot:plev,i0:i1,j0:j1) = op (kbot:plev,i0:i1,j0:j1) ! opnm_out(kbot:plev,i0:i1,j0:j1) = opnm(kbot:plev,i0:i1,j0:j1) ! return ! !------------------------ End third latitude scan --------------------- ! kprint = kbot+5 ! call printpoles(explicit(kbot:plev,i0:i1,j0:j1),kprint,kbot,plev,i0,i1,j0,j1,'EXPLICIT' ) ! call printpoles(p_coeff (kbot:plev,i0:i1,j0:j1),kprint,kbot,plev,i0,i1,j0,j1,'P_COEFF' ) ! call printpoles(q_coeff (kbot:plev,i0:i1,j0:j1),kprint,kbot,plev,i0,i1,j0,j1,'Q_COEFF' ) ! call printpoles(r_coeff (kbot:plev,i0:i1,j0:j1),kprint,kbot,plev,i0,i1,j0,j1,'R_COEFF' ) ! ! Set poles for selected diagnostics: ! call setpoles(diag26(kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! EXPLICITa call setpoles(diag27(kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! EXPLICITb call setpoles(diag2 (kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! EXPLICIT0 call setpoles(diag3 (kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! EXPLICIT1 call setpoles(diag6 (kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! DJINT ! ! All tasks have global 2d bmod2. ! bmod2 was set by sub magfield (getapex.F90) ! allocate(bmod2(0:nlonp1,jspole-1:jnpole+1)) ! Copy bmod2 poles to diagnostic array. ! do i=i0,i1 do k=kbot,plev diag25(k,i,j0) = bmod2(i,j0) diag25(k,i,j1) = bmod2(i,j1) enddo enddo ! call savefld_waccm_switch(diag25,'BMOD2' ,plev,i0,i1,j0,j1) ! ! Assign polar values to coefficients for trsolv. ! call setpoles(explicit(kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) call setpoles(p_coeff (kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) call setpoles(q_coeff (kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) call setpoles(r_coeff (kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! ! Call solver, defining O+ output op_out: ! ! Its best not to call this unless the coefficients and explicit terms ! have been properly set in the third latitude scan above (e.g., during ! "goto 300" debugging above, where the coeffs may not have been calculated). ! if (calltrsolv) then do lat=j0,j1 call trsolv(p_coeff (kbot:plev,i0:i1,lat), & q_coeff (kbot:plev,i0:i1,lat), & r_coeff (kbot:plev,i0:i1,lat), & explicit(kbot:plev,i0:i1,lat), & op_out (kbot:plev,i0:i1,lat), & kbot,plev,kbot,plev,i0,i1,lat) enddo ! call savefld_waccm_switch(op_out,'OP_SOLVE' ,plev,i0,i1,j0,j1) else ! trsolv not called (debug only) op_out (kbot:plev,i0:i1,j0:j1) = op (kbot:plev,i0:i1,j0:j1) opnm_out(kbot:plev,i0:i1,j0:j1) = opnm(kbot:plev,i0:i1,j0:j1) endif ! calltrsolv ! ! Write fields from third latitude scan to waccm history: ! ! call savefld_waccm_switch(explicit,'EXPLICIT',plev,i0,i1,j0,j1) ! non-zero at ubc ! call savefld_waccm_switch(p_coeff ,'P_COEFF' ,plev,i0,i1,j0,j1) ! zero at ubc? ! call savefld_waccm_switch(q_coeff ,'Q_COEFF' ,plev,i0,i1,j0,j1) ! non-zero at ubc ! call savefld_waccm_switch(r_coeff ,'R_COEFF' ,plev,i0,i1,j0,j1) ! is set zero at ubc ! call savefld_waccm_switch(bdotdh_diff(:,i0:i1,j0:j1), 'BDOTDH_DIFF',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag1 ,'BDZDVB_OP',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag2 ,'EXPLICIT0',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag26,'EXPLICITa',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag27,'EXPLICITb',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag3 ,'EXPLICIT1',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag4 ,'TPHDZ0' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag5 ,'TPHDZ1' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag6 ,'DJINT' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag7 ,'DIVBZ' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag8 ,'HDZMBZ' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag9 ,'HDZPBZ' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag10,'EXPLICIT2',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag11,'P_COEFF0' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag12,'Q_COEFF0' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag13,'R_COEFF0' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag14,'BDOTU' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag15,'P_COEFF1' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag16,'Q_COEFF1' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag17,'R_COEFF1' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag18,'EXPLICIT3',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag19,'P_COEFF2' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag20,'Q_COEFF2' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag21,'R_COEFF2' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag22,'P_COEFF0a',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag23,'Q_COEFF0a',plev,i0,i1,j0,j1) ! call savefld_waccm_switch(diag24,'R_COEFF0a',plev,i0,i1,j0,j1) ! ! For debug: return before fourth latitude scan. ! Be careful here: opnm_out has not yet been defined. ! write(iulog,"('oplus returning before fourth lat scan')") ! return !------------------------------------------------------------------------ ! ! Filter O+ output from solver: ! (TIMEGCM calls both filters, whereas TIEGCM calls only filter2) ! ! call filter1_op(op_out(kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! call filter2_op(op_out(kbot:plev,i0:i1,j0:j1),kbot,plev,i0,i1,j0,j1) ! !----------------------- Begin fourth latitude scan --------------------- ! do lat=j0,j1 do i=i0,i1 do k=kbot,plev opnm_out(k,i,lat) = dtsmooth*op(k,i,lat)+dtsmooth_div2* & (opnm(k,i,lat)+op_out(k,i,lat)) enddo enddo ! ! Insure non-negative O+ output: do i=i0,i1 do k=kbot,plev if (op_out (k,i,lat) < 1.e-5_r8) op_out (k,i,lat) = 1.e-5_r8 if (opnm_out(k,i,lat) < 1.e-5_r8) opnm_out(k,i,lat) = 1.e-5_r8 enddo ! k=lev0,lev1-1 enddo ! i=lon0,lon1 enddo ! lat=lat0,lat1 ! ! Save O+ output to WACCM history (cm^3): ! call savefld_waccm_switch(op_out (:,i0:i1,j0:j1),'OP_OUT' ,plev,i0,i1,j0,j1) ! call savefld_waccm_switch(opnm_out(:,i0:i1,j0:j1),'OPNM_OUT',plev,i0,i1,j0,j1) end subroutine oplus_xport !----------------------------------------------------------------------- subroutine oplus_flux(opflux,lon0,lon1,lat0,lat1) ! ! Calculate O+ number flux for sub oplus_xport. ! Flux is returned in opflux(lon0:lon1,lat0:lat1). ! ! alatm: geomagnetic latitude at each geographic grid point (radians) use getapex,only: alatm ! (nlonp1,jspole:jnpole) ! ! Args: integer,intent(in) :: lon0,lon1,lat0,lat1 real(r8),intent(out) :: opflux(lon0:lon1,lat0:lat1) ! ! Local: integer :: i,j,lat real(r8),dimension(lon0:lon1,lat0:lat1) :: chi ! solar zenith angle real(r8),parameter :: & phid = 2.0e8_r8, & phin = -2.0e8_r8, & ! phin = 0._r8, & ppolar = 0._r8 real(r8) :: a(lon0:lon1) real(r8) :: fed(lon0:lon1) real(r8) :: fen(lon0:lon1) ! ! Set some paramaters: pi = 4._r8*atan(1._r8) rtd = 45._r8/atan(1._r8) ! ! Sub get_zenith calls sub zenith (..cam/src/physics/cam/zenith.F90) call get_zenith(chi,lon0,lon1,lat0,lat1) ! ! Latitude scan: do j=lat0,lat1 ! ! Longitude loop: do i=lon0,lon1 if (abs(alatm(i,j))-pi/24._r8>=0._r8) then a(i) = 1._r8 else a(i)=.5_r8*(1._r8+sin(pi*(abs(alatm(i,j))-pi/48._r8)/(pi/24._r8))) if (a(i) < 0.05_r8) a(i) = 0.05_r8 endif fed(i) = phid*a(i) fen(i) = phin*a(i) if (chi(i,j)-0.5_r8*pi >= 0._r8) then opflux(i,j) = fen(i) else opflux(i,j) = fed(i) endif if ((chi(i,j)*rtd-80._r8)*(chi(i,j)*rtd-100._r8) < 0._r8) then opflux(i,j) = .5_r8*(fed(i)+fen(i))+.5_r8*(fed(i)-fen(i))* & cos(pi*(chi(i,j)*rtd-80._r8)/20._r8) endif ! ! Add ppolar if magnetic latitude >= 60 degrees: ! QUESTION: is the 60 deg here related to critical angles crit(2) in tiegcm? ! 3/15/15: opflux is comparable to tiegcm. ! if (abs(alatm(i,j))-pi/3._r8 >= 0._r8) & opflux(i,j) = opflux(i,j)+ppolar enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 ! end subroutine oplus_flux !----------------------------------------------------------------------- subroutine get_zenith(chi,i0,i1,j0,j1) ! ! Get solar zenith angle chi(i0:i1,j0:j1) (radians) ! Subroutine zenith returns cos(chi) at each (i,j) ! Note glon(i0:i1) from edyn_init is in -180 -> +180 (TIEGCM mode) ! use time_manager ,only: get_curr_calday use edyn_geogrid,only : glon,glat use orbit, only : zenith ! ! Args: integer,intent(in) :: i0,i1,j0,j1 real(r8),intent(out) :: chi(i0:i1,j0:j1) ! ! Local: integer :: i,j real(r8) :: dtr,calday real(r8) :: cosZenAngR(1) dtr = pi/180._r8 calday = get_curr_calday() ! fractional day of year do j=j0,j1 do i=i0,i1 call zenith(calday,dtr*glat(j),dtr*glon(i),cosZenAngR,1) chi(i,j) = acos(cosZenAngR(1)) enddo ! write(iulog,"('get_zenith: j=',i4,' glat(j)=',f8.3,' calday=',i4,' chi(:,j)=',/,(6es12.4))") & ! j,glat(j),calday,chi(:,j) enddo end subroutine get_zenith !----------------------------------------------------------------------- subroutine divb(dvb,i0,i1,j0,j1) ! ! Evaluate divergence of B, the unit magnetic field vector. ! (all processors have the full global 2d field) ! ! Args: integer,intent(in) :: i0,i1,j0,j1 real(r8),intent(out) :: dvb(i0:i1,j0:j1) ! ! Local: integer :: i,j,jm1,jp1 real(r8),parameter :: re = 6.37122e8 ! earth radius (cm) dvb = 0._r8 call savefld_waccm(bx(i0:i1,j0:j1),'OPLUS_BX',1,i0,i1,j0,j1) call savefld_waccm(by(i0:i1,j0:j1),'OPLUS_BY',1,i0,i1,j0,j1) call savefld_waccm(bz(i0:i1,j0:j1),'OPLUS_BZ',1,i0,i1,j0,j1) call savefld_waccm(bmod2(i0:i1,j0:j1),'OPLUS_BMAG',1,i0,i1,j0,j1) ! ! Note re is in cm. ! (bx,by,bz are set by sub magfield (getapex.F90)) ! (dphi,dlamda, and cs are set by sub set_geogrid (edyn_init.F90)) ! do j=j0,j1 jm1 = j-1 jp1 = j+1 do i=i0,i1 dvb(i,j) = (((bx(i+1,j)-bx(i-1,j))/(2._r8*dlamda)+ & (cs(jp1)*by(i,jp1)-cs(jm1)*by(i,jm1))/(2._r8*dphi))/ & cs(j)+2._r8*bz(i,j))/re enddo ! i=i0,i1 enddo ! j=j0,j1 end subroutine divb !----------------------------------------------------------------------- subroutine rrk(t,rms,ps1,ps2,n2,tr,ans,lon0,lon1,lev0,lev1,lat) ! ! Returns ambipolar diffusion coefficient in ans. ! ! Args: integer,intent(in) :: lon0,lon1,lev0,lev1,lat real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: & t,rms,ps1,ps2,n2,tr real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: ans ! ! Local: ! integer :: k,i ! do i=lon0,lon1 do k=lev0,lev1-1 ans(k,i) = 1.42E17_r8*boltz*t(k,i)/(p0*expz(k)*.5_r8*(rms(k,i)+ & rms(k+1,i))*(ps2(k,i)*rmassinv_o1*sqrt(tr(k,i))*(1._r8-0.064_r8* & log10(tr(k,i)))**2*colfac+18.6_r8*n2(k,i)*rmassinv_n2+18.1_r8* & ps1(k,i)*rmassinv_o2)) enddo ! k=lev0,lev1 ans(lev1,i) = ans(lev1-1,i) ! should not need to do this enddo ! i=lon0,lon1 ! ! Cap ambipolar diffusion coefficient in ans. ! where( ans(:,:) > 1.5e+8_r8 ) ans(:,:) = 1.5e+8_r8 endwhere end subroutine rrk !----------------------------------------------------------------------- integer function find_kbot(zg,k0,k1,i0,i1,j0,j1,zbot) use edyn_mpi,only: ntask,mytid ! ! Locate k index to pressure level closest to height zg, ! choose the min, and broadcast the index to all tasks. ! #include ! ! Args: integer,intent(in) :: k0,k1,i0,i1,j0,j1 real(r8),intent(in) :: & zg(k0:k1,i0:i1,j0:j1), & ! geopotential height (cm) zbot ! target height (km) ! ! Local: integer :: k,i,j,ktmp real(r8) :: fmin0,fmax0,fmin1,fmax1 integer :: ktmps(ntask),ier find_kbot = 0 ktmp = 0 ! do k=k0,k1-1 fmin0 = minval(zg(k,:,:) *1.e-5_r8) ! minimum km at level k in this subdomain fmax0 = maxval(zg(k,:,:) *1.e-5_r8) ! maximum km at level k in this subdomain fmin1 = minval(zg(k+1,:,:)*1.e-5_r8) ! minimum km at level k+1 in this subdomain fmax1 = maxval(zg(k+1,:,:)*1.e-5_r8) ! maximum km at level k+1 in this subdomain if (fmin0 < zbot.and.fmax0 > zbot) then ktmp = k exit endif if ((fmax0 > zbot.and.fmin1 < zbot).or.(fmin1 < zbot.and.fmax1 > zbot)) then ktmp = k+1 exit endif enddo ! k=k0,k1-1 if (ktmp > 0) then write(iulog,"('Set ktmp=',i4,' zbot=',f9.2,' fmin0,fmax0=',2es12.4,' fmin1,fmax1=',2es12.4)") & ktmp,zbot,fmin0,fmax0,fmin1,fmax1 else write(iulog,"('>>> Could not find kbot: zbot=',f9.2,' km')") zbot call endrun('find_kbot in oplus') endif ! ! Assuming all processors should use the same kbot, gather ktmps from ! the subdomains to the root task, take the minimum of all processor ktmps ! (assuming its best to start integrations slightly below zbot rather ! than above), and scatter it back out to the processors. ! call mpi_gather(ktmp,1,MPI_INTEGER,ktmps,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier) if (ier /= 0) then write(iulog,"('>>> find_kbot: error return from mpi_gather ier=',i5)") ier call endrun endif if (mytid==0) then ktmp = minval(ktmps) ktmps(:) = ktmp endif call mpi_scatter(ktmps,1,MPI_INTEGER,ktmp,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier) if (ier /= 0) then write(iulog,"('>>> find_kbot: error return from mpi_scatter ier=',i5)") ier call endrun endif find_kbot = ktmp ! final index to be used by all processors end function find_kbot !----------------------------------------------------------------------- subroutine diffus(tp,en,hj,ans,i0,i1,lev0,lev1,lat) ! kbot,plev ! Evaluates ans = (d/(h*dz)*tp+m*g/r)*en ! Remember: "bot2top": lev0=kbot=bottom, lev1=plev=top ! ! Args: integer :: i0,i1,lev0,lev1,lat real(r8),dimension(lev0:lev1,i0:i1),intent(in) :: tp,en,hj real(r8),dimension(lev0:lev1,i0:i1),intent(out) :: ans ! ! Local: integer :: i,k real(r8) :: mgr mgr = rmass_op*grav_cm/gask do i=i0,i1 do k=lev0,lev1-2 ans(k+1,i) = 1._r8/(2._r8*hj(k+1,i)*dzp)*(tp(k+2,i)*en(k+2,i)- & tp(k,i)*en(k,i))+mgr*en(k+1,i) enddo ! write(iulog,"('diffus: lat=',i4,' i=',i4,' ans(lev0:lev1-1,i)=',2es12.4)") & ! lat,i,minval(ans(lev0:lev1-1,i)),maxval(ans(lev0:lev1-1,i)) enddo ! ! Upper and lower boundaries: ! do i=i0,i1 ! ! Upper boundary: ans(lev1,i) = 1._r8/(hj(lev1,i)*dzp)*(tp(lev1,i)*en(lev1,i)- & tp(lev1-1,i)*en(lev1-1,i))+mgr*en(lev1,i) ! ! Lower boundary: ans(lev0,i) = 1._r8/(hj(lev0,i)*dzp)*(tp(lev0+1,i)*en(lev0+1,i)- & tp(lev0,i)*en(lev0,i))+mgr*en(lev0,i) enddo end subroutine diffus !----------------------------------------------------------------------- subroutine bdotdh(phijm1,phij,phijp1,ans,lon0,lon1,lev0,lev1,lat) ! ! Evaluates ans = (b(h)*del(h))*phi ! ! Args: integer,intent(in) :: lon0,lon1,lev0,lev1,lat real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: phijm1,phijp1 real(r8),dimension(lev0:lev1,lon0-2:lon1+2),intent(inout) :: phij ! why intent(inout)? real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: ans ! ! Local: integer :: k,i ! ! Note phij longitude dimension is lon0-2:lon1+2 (only i-1 and i+1 are used). ! Halo longitudes i-1 and i+1 must have been set before this routine is ! called. ('by' is use-associated above) ! do i=lon0,lon1 do k=lev0,lev1 ans(k,i) = 1._r8/re*(bx(i,lat)/(cs(lat)*2._r8*dlamda)* & (phij(k,i+1)-phij(k,i-1))+by(i,lat)* & (phijp1(k,i)-phijm1(k,i))/(2._r8*dphi)) enddo ! k=lev0,lev1 enddo ! i=lon0,lon1 ! end subroutine bdotdh !----------------------------------------------------------------------- subroutine bdzdvb(phi,dvb,h,ans,lev0,lev1,lon0,lon1,lat) ! ! Evaluates ans = (bz*d/(h*dz)+divb)*phi ! ! Args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat real(r8),intent(in) :: dvb(lon0:lon1) real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: phi,h real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: ans ! ! Local: integer :: k,i ! do i=lon0,lon1 do k=lev0+1,lev1-1 ans(k,i) = bz(i,lat)/(2._r8*h(k,i)*dzp)*(phi(k+1,i)-phi(k-1,i))+ & dvb(i)*phi(k,i) enddo ! k=lev0+1,lev1-1 enddo ! i=lon0,lon1 ! ! Upper and lower boundaries: do i=lon0,lon1 ans(lev1,i) = bz(i,lat)/(h(lev1,i)*dzp)*(phi(lev1,i)- & phi(lev1-1,i))+dvb(i)*phi(lev1,i) ans(lev0,i) = bz(i,lat)/(h(lev0,i)*dzp)* & (phi(lev0+1,i)-phi(lev0,i))+dvb(i)*phi(lev0,i) enddo ! i=lon0,lon1 end subroutine bdzdvb !----------------------------------------------------------------------- subroutine trsolv(a,b,c,f,x,lev0,lev1,k1,k2,lon0,lon1,lat) ! ! Tri-diagonal solver. ! a(k,i)*x(k-1,i) + b(k,i)*x(k,i) + c(k,i)*x(k+1,i) = f(k,i) ! implicit none ! ! Args: integer,intent(in) :: lev0,lev1,k1,k2,lon0,lon1,lat real(r8),dimension(lev0:lev1,lon0:lon1),intent(in) :: & a, & ! input coefficients b, & ! input coefficients c, & ! input coefficients f ! input RHS real(r8),dimension(lev0:lev1,lon0:lon1),intent(out) :: & x ! output ! ! Local: integer :: k,kk,i real(r8),dimension(lev0:lev1,lon0:lon1) :: & w1,w2,w3 ! work arrays ! ! Lower boundary (W(K1)=B(K1): do i=lon0,lon1 w1(lev0,i) = b(lev0,i) enddo ! ! Set up work arrays: do i=lon0,lon1 do k=k1+1,k2 ! ! W(KF+K-1)=C(K-1)/W(K-1): w2(k-1,i) = c(k-1,i) / w1(k-1,i) ! ! W(K)=A(K)*W(KF+K-1) w1(k,i) = a(k,i) * w2(k-1,i) ! ! W(K)=B(K)-W(K) w1(k,i) = b(k,i) - w1(k,i) enddo ! k=k1+1,k2 enddo ! i=lon0,lon1 ! ! Lower boundary (W(2*KF+K1)=F(K1)/W(K1)): do i=lon0,lon1 w3(k1,i) = f(k1,i) / w1(k1,i) enddo ! do i=lon0,lon1 do k=k1+1,k2 ! ! W(2*KF+K)=A(K)*W(2*KF+K-1) w3(k,i) = a(k,i) * w3(k-1,i) ! ! W(2*KF+K)=F(K)-W(2*KF+K) w3(k,i) = f(k,i) - w3(k,i) ! ! W(2*KF+K)=W(2*KF+K)/W(K) w3(k,i) = w3(k,i) / w1(k,i) enddo ! k=k1+1,k2 enddo ! i=lon0,lon1 ! ! Upper boundary (X(K2)=W(2*KF+K2)): do i=lon0,lon1 x(k2,i) = w3(k2,i) enddo ! ! Back substitution: do i=lon0,lon1 do kk=k1+1,k2 k = k1+k2-kk ! k2-1,k1,-1 ! ! X(K)=W(KF+K)*X(K+1) x(k,i) = w2(k,i) * x(k+1,i) ! ! X(K)=W(2*KF+K)-X(K) x(k,i) = w3(k,i) - x(k,i) enddo ! k=k1+1,k2 enddo end subroutine trsolv !----------------------------------------------------------------------- subroutine printpoles(f,klev,k0,k1,i0,i1,j0,j1,name) ! ! call printpoles(bdotdh_op(kbot:plev,i0:i1,j0:j1),kprint,kbot,plev,i0,i1,j0,j1,'BDOTDH_OP' ) ! ! Args: integer,intent(in) :: klev,k0,k1,i0,i1,j0,j1 real(r8),intent(in) :: f(k0:k1,i0:i1,j0:j1) character(len=*),intent(in) :: name ! ! Print values at the poles at klev: if (j0==1) then write(iulog,"(/,'printpoles ',a,' spole: klev=',i4,' f(klev,i0:i1,j0)=',/,(8es12.4))") & name,klev,f(klev,i0:i1,j0) endif if (j1==plat) then write(iulog,"(/,'printpoles ',a,' npole: klev=',i4,' f(klev,i0:i1,j1)=',/,(8es12.4))") & name,klev,f(klev,i0:i1,j1) endif end subroutine printpoles !----------------------------------------------------------------------- subroutine filter1_op(f,k0,k1,i0,i1,j0,j1) ! ! Polar fft filter, option 1 (see filter.F90). ! use filter_module,only: filter1,filter2,kut1 use edyn_mpi ,only: mp_gatherlons_f3d,mytidi use edyn_mpi ,only: mp_scatterlons_f3d ! ! Args: integer,intent(in) :: k0,k1,i0,i1,j0,j1 real(r8),intent(inout) :: f(k0:k1,i0:i1,j0:j1) ! ! Local: integer :: i,j,k,nlevs,kprint real(r8) :: fik(plon,k1-k0+1) type(array_ptr_type) :: fkij(1) ! fkij(1)%ptr(k1-k0+1,plon,j0:j1) nlevs = k1-k0+1 ! ! Define lons in fkij from current task subdomain: ! allocate(fkij(1)%ptr(nlevs,plon,j0:j1)) do j=j0,j1 do i=i0,i1 do k=k0,k1 ! kbot,plev fkij(1)%ptr(k-k0+1,i,j) = f(k,i,j) enddo enddo enddo ! ! Gather longitudes into tasks in first longitude column of task table ! (leftmost of each j-row) for global fft. (i.e., tasks with mytidi==0 ! gather lons from other tasks in that row). This includes all latitudes. ! call mp_gatherlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1) ! ! Only leftmost tasks at each j-row of tasks does the global filtering: ! if (mytidi==0) then ! ! Define 2d array with all longitudes for filter at each latitude: ! latscan: do j=j0,j1 if (kut1(j) >= plon/2) cycle latscan do i=1,plon do k=k0,k1 fik(i,k-k0+1) = fkij(1)%ptr(k-k0+1,i,j) enddo enddo ! ! Remove wave numbers > kut(lat): ! call filter1(fik,1,nlevs,kut1(j),j) ! ! Return filtered array to fkij: ! do i=1,plon do k=k0,k1 fkij(1)%ptr(k-k0+1,i,j) = fik(i,k-k0+1) enddo enddo ! i=1,plon enddo latscan ! j=j0,j1 endif ! mytidi==0 ! ! Now leftmost task at each j-row must redistribute filtered data ! back to other tasks in the j-row (mytidi>0,mytidj) (includes latitude): ! call mp_scatterlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1) ! ! Return filtered array to inout field at task subdomain: ! do j=j0,j1 do i=i0,i1 do k=k0,k1 f(k,i,j) = fkij(1)%ptr(k-k0+1,i,j) enddo enddo enddo deallocate(fkij(1)%ptr) end subroutine filter1_op !----------------------------------------------------------------------- subroutine filter2_op(f,k0,k1,i0,i1,j0,j1) use filter_module,only: filter1,filter2,kut2 use edyn_mpi ,only: mp_gatherlons_f3d,mytidi use edyn_mpi ,only: mp_scatterlons_f3d ! ! Args: integer,intent(in) :: k0,k1,i0,i1,j0,j1 real(r8),intent(inout) :: f(k0:k1,i0:i1,j0:j1) ! ! Local: integer :: i,j,k,kprint,nlevs real(r8) :: fik(plon,k1-k0+1) type(array_ptr_type) :: fkij(1) ! fkij(1)%ptr(k1-k0+1,plon,j0:j1) nlevs = k1-k0+1 ! ! Define lons in fkij from current task subdomain: ! allocate(fkij(1)%ptr(nlevs,plon,j0:j1)) do j=j0,j1 do i=i0,i1 do k=k0,k1 fkij(1)%ptr(k-k0+1,i,j) = f(k,i,j) enddo enddo enddo ! ! Gather longitudes into tasks in first longitude column of task table ! (leftmost of each j-row) for global fft. (i.e., tasks with mytidi==0 ! gather lons from other tasks in that row). This includes all latitudes. ! call mp_gatherlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1) ! ! Only leftmost tasks at each j-row of tasks does the global filtering: ! if (mytidi==0) then ! ! Define 2d array with all longitudes for filter at each latitude: ! do j=j0,j1 do i=1,plon do k=k0,k1 fik(i,k-k0+1) = fkij(1)%ptr(k-k0+1,i,j) enddo enddo ! ! Remove wave numbers > kut(lat): ! call filter2(fik,1,nlevs,kut2(j),j) ! ! Return filtered array to fkij: ! do i=1,plon do k=k0,k1 fkij(1)%ptr(k-k0+1,i,j) = fik(i,k-k0+1) enddo enddo ! i=1,plon enddo ! j=j0,j1 endif ! mytidi==0 ! ! Now leftmost task at each j-row must redistribute filtered data ! back to other tasks in the j-row (mytidi>0,mytidj) (includes latitude): ! call mp_scatterlons_f3d(fkij,1,nlevs,i0,i1,j0,j1,1) ! ! Return filtered array to inout field at task subdomain: do j=j0,j1 do i=i0,i1 do k=k0,k1 f(k,i,j) = fkij(1)%ptr(k-k0+1,i,j) enddo enddo enddo deallocate(fkij(1)%ptr) end subroutine filter2_op !----------------------------------------------------------------------- end module oplus