module pdynamo_module ! ! This software is part of the NCAR TIE-GCM. Use is governed by the ! Open Source Academic Research License Agreement contained in the file ! tiegcmlicense.txt. ! ! Parallel dynamo module. ! ! use params_module,only: nmlevp1,nmlon,nmlonp1,nlon,nlonp4,nlat, | nmlat,gmlat,nmlath,nlev,nlevp1,mlev0,mlev1,nmlev_diff,dlat,spval use addfld_module,only: addfld use mpitime_module,only: mpi_timer use mpi_module,only: mlat0,mlat1,mlon0,mlon1,lon0,lon1,lat0,lat1, | mp_mag_periodic_f2d use heelis_module,only: phihm use init_module,only: istep use input_module,only: current_pg,current_kq implicit none ! ! 3d pointers to fields regridded to magnetic subdomains (i,j,k): ! (mlon0:mlon1,mlat0:mlat1,mlev0:mlev1) real,pointer,dimension(:,:,:),save :: ped_mag,hal_mag,zpot_mag, | adotv1_mag, adotv2_mag, je1pg_mag, je2pg_mag ! ! 2d pointers to fields on magnetic subdomains (i,j): ! (mlon0:mlon1,mlat0:mlat1) real,pointer,dimension(:,:),save :: sini_mag,adota1_mag, | adota2_mag,a1dta2_mag,be3_mag ! ! 2d coefficients and RHS terms for PDE on magnetic subdomains ! including halo points: real,allocatable,dimension(:,:),save :: | zigm11, ! sigma11*cos(theta0) | zigmc, ! sigmac | zigm2, ! sigma2 | zigm22, ! sigma22/cos(theta0) | rim1,rim2, | rhs, ! right-hand side | phimsolv, ! solution direct from solver (nhem only) | phim2d ! solution with phihm and both nhem and shem #if defined(INTERCOMM) || defined(CISMAH) real,allocatable,dimension(:,:),save :: | zigm1, nsrhs #endif ! ! 3d potential and electric field on mag subdomains (see sub pthreed): ! (mlon0:mlon1,mlat0:mlat1,mlev0:mlev1) real,allocatable,dimension(:,:,:),save :: | phim3d, ! 3d electric potential | ed13d,ed23d, ! 3d electric field for current calculations | ephi3d, ! 3d eastward electric field | elam3d, ! 3d equatorward electric field | emz3d, ! 3d upward electric field | zpotm3d ! 3d geopotential (values at all levels) ! ! 3d electric field on geographic subdomains (see sub pefield): ! (nlevp1,lon0-2,lon1+2,lat0:lat1) real,allocatable,dimension(:,:,:),save :: ex,ey,ez ! ! 3d electric potential on geographic subdomains (lon0:lon1,lat0:lat1,nlevp1) ! This will be regridded from phim3d for output to history files. real,allocatable,dimension(:,:,:),save :: phig3d ! ! coef_rhs and coef_cofum are allocated without halo points: real,allocatable,dimension(:,:),save :: | coef_rhs ! c0(:,:,10) at subdomains (10th only, i.e., rhs) real,allocatable,dimension(:,:,:),save :: coef_cofum ! cofum(mlon0:1,mlat0:1,1-9) ! ! Fields at mag equator (mlon0:mlon1,mlev0:mlev1): real,allocatable,dimension(:,:),save :: | ped_meq, hal_meq, adotv1_meq, adotv2_meq, ! (mlon0:mlon1,mlev0:mlev1) | je1pg_meq, je2pg_meq real,allocatable,dimension(:,:,:),save :: fmeq_out ! (mlon0:mlon1,mlev0:mlev1,6) real,allocatable,dimension(:,:,:,:),save :: fmeq_in ! (mlon0:mlon1,mlat0:mlat1,mlev0:mlev1,6) integer :: jlateq ! latitude index of mag equator real,allocatable,dimension(:,:) :: ressolv ! mag subdomains w/o halos real :: l2norm ! ! Global longitude values near mag equator and poles for complete_integrals and rhs. ! The nf2d 6 fields are: zigm11,zigm22,zigmc,zigm2,rim1,rim2 #if defined(INTERCOMM) || defined(CISMAH) integer,parameter :: nf2d=7 ! 7 2d fields for CISM #else integer,parameter :: nf2d=6 ! 6 2d fields #endif real :: feq_jpm1(nmlonp1,2,nf2d) ! 6 fields at 2 lats (eq-1, eq+1) real :: fpole_jpm2(nmlonp1,4,nf2d) ! fields at S pole+1,2 and N pole-1,2 ! ! Declarations for remaining serial part of dynamo: ! (see sub stencils, solver, and related subroutines) ! ! Global 2d fields for root task to complete non-parallel part of dynamo. ! Gathered from subdomains by in sub gather_pdyn. real,dimension(nmlonp1,nmlat) :: phim_glb, ! phim_glb is used by noso_crrt | zigm11_glb,zigm22_glb,zigmc_glb,zigm2_glb,rhs_glb real,dimension(nmlonp1,nmlat,2) :: rim_glb ! pde solver output real,dimension(0:nmlonp1,0:nmlat+1) :: phisolv ! ! Dimensions of the 5 grid resolutions for the multi-grid PDE: integer,parameter :: | nmlon0=nmlon+1, | nmlat0=(nmlat+1)/2, | nmlon1=(nmlon0+1)/2, | nmlat1=(nmlat0+1)/2, | nmlon2=(nmlon1+1)/2, | nmlat2=(nmlat1+1)/2, | nmlon3=(nmlon2+1)/2, | nmlat3=(nmlat2+1)/2, | nmlon4=(nmlon3+1)/2, | nmlat4=(nmlat3+1)/2 ! ! Unmodified coefficients for using modified mudpack: real,dimension(nmlon0,nmlat0,9) :: cofum ! ! Space needed for descretized coefficients of of dynamo pde at all ! 5 levels of resolution: ! integer,parameter :: | ncee=10*nmlon0*nmlat0+9*(nmlon1*nmlat1+nmlon2*nmlat2+nmlon3* | nmlat3+nmlon4*nmlat4) ! ! Coefficients are stored in 1-d array cee(ncee) ! cee transmits descretized dynamo PDE coefficients to the multi-grid ! mudpack solver. (cee was formerly in ceee.h) ! The common block /cee_com/ is retained from earlier versions because ! of the equivalencing below of coefficient arrays c0, c1, etc. ! real :: cee(ncee) common/cee_com/ cee ! ! The following parameters nc0,nc1,... are pointers to the beginning of ! the coefficients for each level of resolution. ! integer,parameter :: | nc0=1, | nc1=nc0+10*nmlon0*nmlat0, | nc2=nc1+9 *nmlon1*nmlat1, | nc3=nc2+9 *nmlon2*nmlat2, | nc4=nc3+9 *nmlon3*nmlat3 ! ! nc(1:6) are pointers to beginning of coefficient blocks at each of ! 5 levels of resolution: ! nc(1) = nc0, pointer to coefficients for highest resolution. ! nc(2) = nc1, pointer to coefficients at half the resolution of nc0, ! and so on for nc(3), nc(4), nc(5), etc. ! nc(6) = ncee, the dimension of the entire cee array, containing ! coefficients for all 5 levels of resolution. ! integer :: nc(6) integer :: mlevtop ! for saving mag fields to sechist real :: | c0(nmlon0,nmlat0,10), | c1(nmlon1,nmlat1,9), | c2(nmlon2,nmlat2,9), | c3(nmlon3,nmlat3,9), | c4(nmlon4,nmlat4,9) equivalence | (cee,c0), | (cee(nc1),c1), | (cee(nc2),c2), | (cee(nc3),c3), | (cee(nc4),c4) logical :: debug = .false. real,dimension(nmlonp1,nmlat0) :: pfrac ! NH fraction of potential contains !----------------------------------------------------------------------- subroutine pdynamo use mpi_module,only: mytid,mp_scatter_coeffs,mp_scatter_phim, | mp_mag_halos,mp_gather_pdyn use diags_module,only: mkdiag_PHIM2D ! use solver_module,only: solve_using_pcg ! ! Parallel dynamo. ! Required fields from neutral atmos have been regridded to magnetic ! grid by sub dynamo_inputs (called from advance). ! ! Transform needed fields to geomagnetic coordinates ! Perform field-line integrations ! Evaluate PDE coefficients and RHS ! The PDE is divided by 1/ DT0DTS (in dyncal divided by 1/cos(theta_0) ! Sigma_(phi phi) = zigm11/ rcos0s * dt0dts ! Sigma_(lam lam) = zigm22 * rcos0s / dt0dts ! Sigma_(phi lam) = +-(zigm2-zigmc) ! Sigma_(lam phi) = -+(zigm2+zigmc) ! K_(m phi)^D = rim(1) * dt0dts ! K_(m lam)^D = +-rim(2) * rcos0s ! ! (Sub dynamo_inputs has already been called from advance) ! ! Local: integer :: i,j,n,ier,lonend real :: coefglb_rhs(nmlonp1,nmlat),coefglb_cofum(nmlonp1,nmlat,9) character(len=16) :: cofum_name ! ! Latitude index to mag equator (module data) jlateq = (nmlat+1)/2 ! index to mag equator (module data) ! ! Fieldline integration: call mpi_timer("pdynamo_fieldline_integrals",0,0) ! name,startstop,iprint call fieldline_integrals call mpi_timer("pdynamo_fieldline_integrals",1,0) if (debug) write(6,"('pdynamo after fieldline_integrals')") ! ! Equatorial and polar values, hemisphere folding: call mpi_timer("pdynamo_complete_integrals",0,0) ! name,startstop,iprint call complete_integrals call mpi_timer("pdynamo_complete_integrals",1,0) ! name,startstop,iprint if (debug) write(6,"('pdynamo after complete_integrals')") ! ! Calculate right-hand side on mag subdomains: ! (mag halos are needed in rim1 for rhs calculation) call mp_mag_halos(rim1,mlon0,mlon1,mlat0,mlat1,1,1) call mp_mag_halos(rim2,mlon0,mlon1,mlat0,mlat1,1,1) call mpi_timer("pdynamo_rhspde",0,0) ! name,startstop,iprint call rhspde call mpi_timer("pdynamo_rhspde",1,0) ! name,startstop,iprint if (debug) write(6,"('pdynamo after rhspde')") ! ! Gather needed arrays to root task and complete dynamo serially as in ! original code. Timing is calculated for mp_gather_pdyn, which is ! called by sub gather_pdyn, and is the largest wallclock consumer ! of the entire parallel dynamo code. ! call gather_pdyn if (debug) write(6,"('pdynamo after gather_pdyn')") ! ! Define global mag coeffs arrays, with non-zero in shem indices ! real :: coefglb_rhs(nmlonp1,nmlat),coefglb_cofum(nmlonp1,nmlat,9) coefglb_rhs = 0. ! whole-array init coefglb_cofum = 0. ! whole-array init ! ! Complete dynamo serially by root task: if (mytid==0) then ! root task call mpi_timer("stencils",0,0) call stencils ! set up stencils call mpi_timer("stencils",1,0) if (debug) write(6,"('pdynamo: root task after stencils')") ! call mpi_timer("solver",0,0) call solver(cofum,c0) ! mudpack PDE solver call mpi_timer("solver",1,0) if (debug) write(6,"('pdynamo: root task after solver')") coefglb_rhs (:,1:nmlath) = c0(:,:,10) coefglb_cofum(:,1:nmlath,:) = cofum(:,:,:) endif ! ! rim1 after solver is needed for highlat_poten. rim_glb is distributed ! to subdomains as rim1, and mag halos set. This will overwrite rim1 from ! fieldline_integrals, complete_integrals, etc. ! call mp_scatter_phim(rim_glb(:,:,1),rim1(mlon0:mlon1, | mlat0:mlat1)) if (debug) write(6,"('pdynamo: after mp_scatter_phim')") call mp_mag_halos(rim1,mlon0,mlon1,mlat0,mlat1,1,1) ! ! Add high latitude potential from empirical model (heelis or weimer) ! to solution rim1, defining phim2d on mag subdomains. ! allocate(phim2d(mlon00:mlon11,mlat00:mlat11),stat=istat) ! call mpi_timer("highlat_poten",0,0) call highlat_poten call mpi_timer("highlat_poten",1,0) if (debug) write(6,"('pdynamo: after highlat_poten')") ! write(6,"('pdynamo after highlat_poten: phim2d=',2e12.4)") ! | minval(phim2d(mlon0:mlon1,mlat0:mlat1)), ! | maxval(phim2d(mlon0:mlon1,mlat0:mlat1)) ! ! Gather phim2d to global phim_glb on root task, for use by noso_crrt. ! if (current_kq > 0) then call mp_gather_pdyn(phim2d(mlon0:mlon1,mlat0:mlat1), | mlon0,mlon1,mlat0,mlat1,phim_glb,nmlonp1,nmlat,1) if (debug) write(6,"('pdynamo after mp_gather_pdyn')") endif call mkdiag_PHIM2D('PHIM2D',phim2d(mlon0:mlon1,mlat0:mlat1), | mlon0,mlon1,mlat0,mlat1) ! ! Expand phim2d to phim3d, first setting mag halos in phim2d from ! hightlat_poten. phim3d will then be the final potential from pdynamo. ! call mp_mag_halos(phim2d,mlon0,mlon1,mlat0,mlat1,1,1) if (debug) write(6, | "('pdynamo after mp_mag_halos, before pthreed')") call mpi_timer("pthreed",0,0) call pthreed call mpi_timer("pthreed",1,0) if (debug) write(6,"('pdynamo after pthreed')") ! ! Distribute 2d global phisolv (mudpack solution) to subdomains ! in phimsolv (exclude phim halos in this call). phimsolv is used ! for residual testing by productAX below. ! ! call mp_scatter_phim(phisolv(1:nmlonp1,1:nmlat), ! | phimsolv(mlon0:mlon1,mlat0:mlat1)) ! ! Scatter coeffs from global to subdomains for productAX: ! call mp_scatter_coeffs(coefglb_rhs,coefglb_cofum, ! | coef_rhs,coef_cofum) ! ! matrix vector product: ! ressolv = 0.0 ! mag subdomains w/o halos ! ! Call productAX with phimsolv on subdomains (distributed phisolv), ! and with halos set: ! call productAX(coef_cofum,phimsolv,mlon0,mlon1, ! | mlat0,mlat1,ressolv) ! call solve_using_pcg(phisolv,coef_cofum, ! | coef_rhs(mlon0:mlon1,mlat0:mlat1),10,1e-6, ! | mlon0,mlon1,mlat0,mlat1,nmlon,nmlat) ! calculate residual on subdomain . ! l2norm = 0. ! lonend = mlon1 ! if (mlon1 == nmlon+1) lonend = mlon1-1 . ! do j = mlat0,mlat1 ! do i = mlon0,lonend ! ressolv(i,j) = coef_rhs(i,j)-ressolv(i,j) ! l2norm = l2norm + ressolv(i,j)*ressolv(i,j) ! enddo ! enddo ! write(6,*) 'L2norm (task ',mytid,' subdomain) = ',l2norm if (debug) write(6,"('pdynamo returning')") end subroutine pdynamo !----------------------------------------------------------------------- subroutine dynamo_inputs(un,vn,wn,zpot,ped,hal, | lev0,lev1,lon0,lon1,lat0,lat1) ! ! Provide needed field inputs to the dynamo by regridding the fields ! from geo to mag. ! use params_module,only: dlev,zpbot_dyn,zibot use cons_module,only: h0,dtr use magfield_module,only: zb,bmod use mpi_module,only: mp_mageq,mp_periodic_f3d,mp_periodic_f2d, | mp_magpole_3d use mpi_module,only: tasks,nmagtaski,mytid ! only for check use my_esmf,only: | geo_src_grid, mag_des_grid, esmf_regrid, | esmf_set2d_geo, esmf_get_2dfield, | esmf_set3d_geo, esmf_get_3dfield use my_esmf,only: ! esmf 3d (i,j,k) fields on geo grid | geo_ped, geo_hal, geo_zpot, geo_adotv1, geo_adotv2, | geo_je1pg, geo_je2pg use my_esmf,only: ! esmf 2d (i,j) fields on geo grid | geo_sini, geo_adota1, geo_adota2, geo_a1dta2, geo_be3 use my_esmf,only: ! esmf 3d (i,j,k) fields on mag grid | mag_ped, mag_hal, mag_zpot, mag_adotv1, mag_adotv2, | mag_je1pg, mag_je2pg, esmf_3dfields use my_esmf,only: ! esmf 2d (i,j) fields on mag grid | mag_sini, mag_adota1, mag_adota2, mag_a1dta2, mag_be3 use params_module,only: glat,glon use magpres_g,only: ! (nlevp1,lon0:lon1,lat0:lat1) | je1oD_geo, ! J_e1(pg)/D [A/cm^2] | je2oD_geo ! J_e2(pg)/D [A/cm^2] ! ! Args: integer :: lev0,lev1,lon0,lon1,lat0,lat1,rc,lonbeg,lonend real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(in):: | un, ! neutral zonal velocity (cm/s) | vn, ! neutral meridional velocity (cm/s) | wn ! neutral vertical velocity (cm/s) real,dimension(lev0:lev1,lon0:lon1,lat0:lat1), | intent(inout):: | zpot, ! geopotential height | ped, ! pedersen conductivity (lamdas.F) | hal ! hall conductivity (lamdas.F) ! ! Local: integer :: i,j,k ! ! 3d locals dimensioned mlev0:mlev1 (mlev0 = -2 at 5-deg, -5 at 2.5-deg) real,dimension(mlev0:mlev1,lon0:lon1,lat0:lat1) :: | undyn,vndyn,wndyn,zpotdyn,peddyn,haldyn real,dimension(mlev0:mlev1,lon0:lon1,lat0:lat1) :: | adotv1,adotv2,je1pg,je2pg real,dimension(lon0:lon1,lat0:lat1) :: | sini,adota1,adota2,a1dta2,be3 real,dimension(lon0:lon1,lat0:lat1,4) :: f2d real,parameter :: unitvm(nmlon)=1. ! ! Set constants: ! rl1,rl2 are rates at which sigma1 and sigma2 decay with height ! below bottom of model. ! real,parameter :: | rl1 = 5.e5, | rl2 = 3.e5 ! ! External: real,external :: sddot ! util.F ! ! For esmf_set3d_geo call: real,dimension(mlev0:mlev1,lon0:lon1,lat0:lat1,8) :: f3d character(len=8) :: fnames(8) ! ! write(6,"('Enter dynamo_inputs: un=',2e12.4,' vn=',2e12.4, ! | ' wn=',2e12.4,' zpot=',2e12.4,' ped=',2e12.4,' hal=',2e12.4)") ! | minval(un),maxval(un),minval(vn),maxval(vn), ! | minval(wn),maxval(wn),minval(zpot),maxval(zpot), ! | minval(ped),maxval(ped),minval(hal),maxval(hal) ! ! mlevtop is for saving mag fields with addfld, see below: ! (actually mlevtop one level below the top): ! At 5.0-deg: mlev0,mlev1 = -2,29 mag_field(81,97,mlev0:mlev1) ! At 2.5-deg: mlev0,mlev1 = -5,57 mag_field(81,97,mlev0:mlev1) ! mlevtop = mlev1+2 ! 5.0-deg (mlev1=57) if (dlat == 2.5) mlevtop = mlev1+5 ! 2.5-deg (mlev1=63) ! ! First, transfer lev0:lev1 (and all geo i,j) from inputs to xxxdyn: do j=lat0,lat1 do i=lon0,lon1 do k=lev0,lev1 undyn(k,i,j) = un(k,i,j) vndyn(k,i,j) = vn(k,i,j) wndyn(k,i,j) = wn(k,i,j) zpotdyn(k,i,j) = zpot(k,i,j) peddyn(k,i,j) = ped(k,i,j) haldyn(k,i,j) = hal(k,i,j) je1pg(k,i,j) = je1oD_geo(k,i,j) je2pg(k,i,j) = je2oD_geo(k,i,j) enddo ! k=lev0,lev1 ! ! Extend neutral atmosphere inputs down to mlev0 (zpbot_dyn=-8.5) ! At 5-deg res, this will be k=0,-2, and at 2.5-deg res, k=0,-5. ! Set three equally spaced levels for Z, take U, V, and W ! to be constant, and extrapolate sigmas exponentially. ! ! 3/3/16 btf (as per Astrid and Art): Change extrapolation of Pedersen and Hall: ! from (float(nmlev_diff-1)*rl1)) to (2.*rl1)) for peddyn, and ! from (float(nmlev_diff-1)*rl2)) to (2.*rl2)) for haldyn. ! According to Astrid This will lead to smaller conductivities between ! zp -8.5 to zp -7. ! do k=lev0-1,mlev0,-1 zpotdyn(k,i,j) = h0+float(k+nmlev_diff-1)* | (zpotdyn(lev0,i,j)-h0)/float(nmlev_diff) peddyn(k,i,j) = peddyn(lev0,i,j)*exp((zpotdyn(k,i,j)+ | zpotdyn(k+1,i,j)-zpotdyn(lev0,i,j)-zpotdyn(lev0+1,i,j))/ | (2.*rl1)) haldyn(k,i,j) = haldyn(1,i,j)*exp((zpotdyn(k,i,j)+ | zpotdyn(k+1,i,j)-zpotdyn(lev0,i,j)-zpotdyn(lev0+1,i,j))/ | (2.*rl2)) undyn(k,i,j) = undyn(lev0,i,j) vndyn(k,i,j) = vndyn(lev0,i,j) wndyn(k,i,j) = wndyn(lev0,i,j) je1pg(k,i,j) = 0. je2pg(k,i,j) = 0. enddo enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 ! do k=lev0-1,mlev0,-1 ! write(6,"('dynamo_inputs: k=',i4,' mlev0=',i4,' undyn=',2e12.4, ! | ' vndyn=',2e12.4,' wndyn=',2e12.4,' zpotdyn=',2e12.4, ! | ' peddyn=',2e12.4,' haldyn=',2e12.4)") k,mlev0, ! | minval(undyn(k,:,:)),maxval(undyn(k,:,:)), ! | minval(vndyn(k,:,:)),maxval(vndyn(k,:,:)), ! | minval(wndyn(k,:,:)),maxval(wndyn(k,:,:)), ! | minval(zpotdyn(k,:,:)),maxval(zpotdyn(k,:,:)), ! | minval(peddyn(k,:,:)),maxval(peddyn(k,:,:)), ! | minval(haldyn(k,:,:)),maxval(haldyn(k,:,:)) ! enddo ! ! real,dimension(mlev0:mlev1,lon0:lon1,lat0:lat1) :: do j=lat0,lat1 call addfld('Z_DYN',' ',' ',zpotdyn(1:nlev,:,j), | 'ilev',lev0,lev1-1,'lon',lon0,lon1,j) call addfld('U_DYN',' ',' ',undyn(1:nlev,:,j), | 'lev',lev0,lev1-1,'lon',lon0,lon1,j) call addfld('V_DYN',' ',' ',vndyn(1:nlev,:,j), | 'lev',lev0,lev1-1,'lon',lon0,lon1,j) call addfld('W_DYN',' ',' ',wndyn(1:nlev,:,j), | 'lev',lev0,lev1-1,'lon',lon0,lon1,j) call addfld('PED_DYN',' ',' ',peddyn(1:nlev,:,j), | 'lev',lev0,lev1-1,'lon',lon0,lon1,j) call addfld('HAL_DYN',' ',' ',haldyn(1:nlev,:,j), | 'lev',lev0,lev1-1,'lon',lon0,lon1,j) enddo ! ! Calculate some 2d and 3d fields: call calc_adotv(zpotdyn,undyn,vndyn,wndyn, | adotv1,adotv2,adota1,adota2, | a1dta2,be3,mlev0,mlev1,lon0,lon1,lat0,lat1) ! ! Calculate sini sin(I_m) (zb and bmod are from magfield) ! real,dimension(nlonp1,0:nlatp1) :: xb,bmod ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 do j=lat0,lat1 do i=lonbeg,lonend sini(i,j) = zb(i-2,j)/bmod(i-2,j) ! sin(I_m) enddo if (lon0==1) sini(1:2,j) = | zb(nlon:nlon+1,j)/bmod(nlon:nlon+1,j) if (lon1==nlonp4) sini(nlon+3:nlonp4,j) = | zb(1:2,j)/bmod(1:2,j) enddo ! ! call addfld('BMOD','2d from apex_subs',' ', ! | bmod(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ZB','2d from apex_subs',' ', ! | zb(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0) ! subroutine mp_periodic_f3d(f,lev0,lev1,lon0,lon1,lat0,lat1,nf) ! real,dimension(lev0:lev1,lon0:lon1,lat0:lat1,4) :: f3d ! f3d(:,:,:,1) = peddyn(:,lon0:lon1,lat0:lat1) f3d(:,:,:,2) = haldyn(:,lon0:lon1,lat0:lat1) f3d(:,:,:,3) = adotv1 f3d(:,:,:,4) = adotv2 call mp_periodic_f3d(f3d,mlev0,mlev1,lon0,lon1,lat0,lat1,4) peddyn(:,lon0:lon1,lat0:lat1) = f3d(:,:,:,1) haldyn(:,lon0:lon1,lat0:lat1) = f3d(:,:,:,2) adotv1 = f3d(:,:,:,3) adotv2 = f3d(:,:,:,4) ! ! Set 3d field values on geographic source grid: ! fnames = (/'PED ','HAL ','ZPOT ','SCHT ', | 'ADOTV1 ','ADOTV2 ','JE1PG ','JE2PG '/) f3d(:,:,:,1) = peddyn f3d(:,:,:,2) = haldyn f3d(:,:,:,3) = zpotdyn f3d(:,:,:,4) = adotv1 f3d(:,:,:,5) = adotv2 f3d(:,:,:,6) = je1pg f3d(:,:,:,7) = je2pg esmf_3dfields(1) = geo_ped esmf_3dfields(2) = geo_hal esmf_3dfields(3) = geo_zpot esmf_3dfields(4) = geo_adotv1 esmf_3dfields(5) = geo_adotv2 esmf_3dfields(6) = geo_je1pg esmf_3dfields(7) = geo_je2pg call esmf_set3d_geo(esmf_3dfields,fnames, | f3d,7,mlev0,mlev1,lon0,lon1,lat0,lat1) geo_ped = esmf_3dfields(1) geo_hal = esmf_3dfields(2) geo_zpot = esmf_3dfields(3) geo_adotv1 = esmf_3dfields(4) geo_adotv2 = esmf_3dfields(5) geo_je1pg = esmf_3dfields(6) geo_je2pg = esmf_3dfields(7) peddyn = f3d(:,:,:,1) haldyn = f3d(:,:,:,2) zpotdyn = f3d(:,:,:,3) adotv1 = f3d(:,:,:,4) adotv2 = f3d(:,:,:,5) je1pg = f3d(:,:,:,6) je2pg = f3d(:,:,:,7) ! ! 2d fields need only be calculated in first timestep: if (istep==1) then f2d(:,:,1) = adota1 f2d(:,:,2) = adota2 f2d(:,:,3) = a1dta2 f2d(:,:,4) = be3 call mp_periodic_f2d(f2d,lon0,lon1,lat0,lat1,4) adota1 = f2d(:,:,1) adota2 = f2d(:,:,2) a1dta2 = f2d(:,:,3) be3 = f2d(:,:,4) ! ! Set 2d field values on geographic grid: ! (esmf fields on source grid exclude periodic points) ! call esmf_set2d_geo(geo_sini, geo_src_grid,'SINI ', | sini, lon0,lon1,lat0,lat1) call esmf_set2d_geo(geo_adota1,geo_src_grid,'ADOTA1 ', | adota1,lon0,lon1,lat0,lat1) call esmf_set2d_geo(geo_adota2,geo_src_grid,'ADOTA2 ', | adota2,lon0,lon1,lat0,lat1) call esmf_set2d_geo(geo_a1dta2,geo_src_grid,'A1DTA2 ', | a1dta2,lon0,lon1,lat0,lat1) call esmf_set2d_geo(geo_be3, geo_src_grid,'BE3 ', | be3, lon0,lon1,lat0,lat1) endif ! ! Save original fields on geographic grid to secondary histories: ! do j=lat0,lat1 call addfld('PED_GEO',' ','S/m',ped(:,:,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) call addfld('HAL_GEO',' ','S/m',hal(:,:,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) call addfld('ZPOT_GEO',' ','cm',zpot(:,:,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) ! ! Beware, adotv1,2, etc have levels mlev0:mlev1, e.g., -2:nlevp1, ! so must specify lev0:lev1 ! call addfld('ADOTV1_GEO',' ',' ',adotv1(lev0:lev1,:,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) call addfld('ADOTV2_GEO',' ',' ',adotv2(lev0:lev1,:,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) call addfld('JE1PG_GEO',' ',' ',je1pg(lev0:lev1,:,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) call addfld('JE2PG_GEO',' ',' ',je2pg(lev0:lev1,:,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) enddo ! call addfld('SINI_GEO',' ',' ',sini, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ADOTA1_GEO',' ',' ',adota1, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ADOTA2_GEO',' ',' ',adota2, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('A1DTA2_GEO',' ',' ',a1dta2, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('BE3_GEO',' ',' ',be3, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! ! Regrid 3d geo fields to mag grid: call esmf_regrid(geo_ped ,mag_ped, 'geo2mag',3) call esmf_regrid(geo_hal ,mag_hal, 'geo2mag',3) call esmf_regrid(geo_zpot ,mag_zpot, 'geo2mag',3) call esmf_regrid(geo_adotv1 ,mag_adotv1, 'geo2mag',3) call esmf_regrid(geo_adotv2 ,mag_adotv2, 'geo2mag',3) call esmf_regrid(geo_je1pg ,mag_je1pg, 'geo2mag',3) call esmf_regrid(geo_je2pg ,mag_je2pg, 'geo2mag',3) ! ! Regrid 2d geo fields to mag grid: if (istep==1) then call esmf_regrid(geo_sini ,mag_sini , 'geo2mag',2) call esmf_regrid(geo_adota1 ,mag_adota1, 'geo2mag',2) call esmf_regrid(geo_adota2 ,mag_adota2, 'geo2mag',2) call esmf_regrid(geo_a1dta2 ,mag_a1dta2, 'geo2mag',2) call esmf_regrid(geo_be3 ,mag_be3 , 'geo2mag',2) endif ! ! Define pdynamo module data pointers to the regridded mag fields. ! First arg of esmf_get_field is input esmf field (my_esmf module), ! second arg is output data pointer (pdynamo module) ! (These destination grid fields have periodic points allocated and set) ! ! Get regridded 3d mag fields: ! subroutine esmf_get_3dfield(input field, output fptr, name) ! First arg is input ESMF field, second is output pointer: ! call esmf_get_3dfield(mag_ped ,ped_mag, "PED ") call esmf_get_3dfield(mag_hal ,hal_mag, "HAL ") call esmf_get_3dfield(mag_zpot ,zpot_mag, "ZPOT ") call esmf_get_3dfield(mag_adotv1,adotv1_mag,"ADOTV1 ") call esmf_get_3dfield(mag_adotv2,adotv2_mag,"ADOTV2 ") call esmf_get_3dfield(mag_je1pg ,je1pg_mag, "JE1PG ") call esmf_get_3dfield(mag_je2pg ,je2pg_mag, "JE1PG ") ! ! Get regridded 2d mag fields: ! First arg is input ESMF field, second is output pointer: ! if (istep==1) then call esmf_get_2dfield(mag_sini ,sini_mag , "SINI ") call esmf_get_2dfield(mag_adota1,adota1_mag, "ADOTA1 ") call esmf_get_2dfield(mag_adota2,adota2_mag, "ADOTA2 ") call esmf_get_2dfield(mag_a1dta2,a1dta2_mag, "A1A2M ") call esmf_get_2dfield(mag_be3 ,be3_mag , "BE3 ") endif ! ! Set mag equator: ! fmeq_in are input fields on 3d mag subdomains. ! ! allocate(fmeq_in(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1,6) ! allocate(fmeq_out(mlon0:mlon1,mlev0:mlev1,6),stat=istat) ! allocate(ped_meq(mlon0:mlon1,mlev0:mlev1),stat=istat) ! fmeq_in(:,:,:,1) = ped_mag(:,:,:) fmeq_in(:,:,:,2) = hal_mag(:,:,:) fmeq_in(:,:,:,3) = adotv1_mag(:,:,:) fmeq_in(:,:,:,4) = adotv2_mag(:,:,:) fmeq_in(:,:,:,5) = je1pg_mag(:,:,:) fmeq_in(:,:,:,6) = je2pg_mag(:,:,:) ! ! Tasks w/ mag equator send eq data(i,k) to other tasks in their tidi: ! call mp_mageq(fmeq_in,fmeq_out,6,mlon0,mlon1,mlat0,mlat1, | mlev0,mlev1) ! ! Output arrays now have mag equator data on longitude subdomain ! and full column (mlon0:mlon1,mlev0:mlev1) ! These will be used in fieldline_integrals. ! ped_meq(:,:) = fmeq_out(:,:,1) hal_meq(:,:) = fmeq_out(:,:,2) adotv1_meq(:,:) = fmeq_out(:,:,3) adotv2_meq(:,:) = fmeq_out(:,:,4) je1pg_meq(:,:) = fmeq_out(:,:,5) je2pg_meq(:,:) = fmeq_out(:,:,6) ! ! Save geopotential on magnetic grid in zpotm3d, then ! limit max zpot_mag to h0 for use in fieldline integrals ! and pthreed. It is not necessary to set poles of zpotm3d ! (as in old dynamo for zpotenm3d), since pthreed does not ! reference the poles of zpotm3d). ! ! allocate(zpotm3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) ! do k=mlev0,mlev1 do j=mlat0,mlat1 do i=mlon0,mlon1 zpotm3d(i,j,k) = zpot_mag(i,j,k) if (zpot_mag(i,j,k) < h0) zpot_mag(i,j,k)=h0 enddo enddo enddo ! ! Save mag fields to secondary histories ! ! ZMAG is a mandatory mag field on secondary histories, so this ! addfld call should not be commented. ! (note cannot use mlev0=-2 in addfld call) ! do j=mlat0,mlat1 call addfld('ZMAG','ZMAG from pdynamo','km', | zpot_mag(:,j,mlev0:mlev1-1)*1.e-5,'mlon',mlon0,mlon1, | 'imlev',1,mlevtop,j) enddo do j=mlat0,mlat1 call addfld('PED_MAG','PED on mag grid','S/m', | ped_mag(:,j,mlev0:mlev1-1),'mlon',mlon0,mlon1, | 'mlev',1,mlevtop,j) call addfld('HAL_MAG','HALL on mag grid','S/m', | hal_mag(:,j,mlev0:mlev1-1),'mlon',mlon0,mlon1, | 'mlev',1,mlevtop,j) call addfld('ZPOT_MAG','ZPOT on mag grid','cm', | zpot_mag(:,j,mlev0:mlev1-1),'mlon',mlon0,mlon1, | 'imlev',1,mlevtop,j) call addfld('ADOTV1_MAG','ADOTV1 on mag grid',' ', | adotv1_mag(:,j,mlev0:mlev1-1),'mlon',mlon0,mlon1, | 'mlev',1,mlevtop,j) call addfld('ADOTV2_MAG','ADOTV2 on mag grid',' ', | adotv2_mag(:,j,mlev0:mlev1-1),'mlon',mlon0,mlon1, | 'mlev',1,mlevtop,j) call addfld('JE1PG_MAG','JE1PG on mag grid',' ', | je1pg_mag(:,j,mlev0:mlev1-1),'mlon',mlon0,mlon1, | 'mlev',1,mlevtop,j) call addfld('JE2PG_MAG','JE2PG on mag grid',' ', | je2pg_mag(:,j,mlev0:mlev1-1),'mlon',mlon0,mlon1, | 'mlev',1,mlevtop,j) enddo ! call addfld('SINI_MAG', 'SINI on mag grid',' ', ! | sini_mag,'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ADOTA1_MAG','ADOTA1 on mag grid',' ', ! | adota1_mag,'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ADOTA2_MAG','ADOTA2 on mag grid',' ', ! | adota2_mag,'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('A1DTA2_MAG','A1DTA2 on mag grid',' ', ! | a1dta2_mag,'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('BE3_MAG' ,'BE3 on mag grid',' ', ! | be3_mag,'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! ! if (current_pg > 0) then ! do j=mlat0,mlat1 ! call addfld('JE1PG', ! | 'JE1PG: Eastward current due to plasma pressure and gravity', ! | 'A/cm^2',je1pg_mag(:,j,1:mlev1),'mlon',mlon0,mlon1,'mlev', ! | 1,mlev1,j) ! call addfld('JE2PG', ! | 'JE2PG: Downward current due to plasma pressure and gravity', ! | 'A/cm^2',je2pg_mag(:,j,1:mlev1),'mlon',mlon0,mlon1,'mlev', ! | 1,mlev1,j) ! enddo ! endif end subroutine dynamo_inputs !----------------------------------------------------------------------- subroutine calc_adotv(z,un,vn,wn,adotv1,adotv2,adota1,adota2, | a1dta2,be3,lev0,lev1,lon0,lon1,lat0,lat1) use cons_module,only: r0,h0 use getapex_module,only: dvec ! (nlonp1,nlat,3,2) use getapex_module,only: dddarr ! (nlonp1,nlat) use getapex_module,only: be3arr ! (nlonp1,nlat) use magfield_module,only: alatm ! (nlonp1,0:nlatp1) ! ! Args: ! (actual level args are mlev0,mlev1, and 3d arrays have these dims) ! integer :: lev0,lev1,lon0,lon1,lat0,lat1 real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(in):: | z, ! geopotential height (cm) | un, ! neutral zonal velocity (cm/s) | vn ! neutral meridional velocity (cm/s) real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(in) :: | wn ! vertical velocity (cm/s) real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(out) :: | adotv1, adotv2 real,dimension(lon0:lon1,lat0:lat1),intent(out) :: | adota1, adota2, a1dta2, be3 ! ! Local: integer :: k,i,j,lonbeg,lonend,im2 real :: r0or,rat,sinalat real :: clm2(lon0:lon1,lat0:lat1) ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 adotv1 = 0. ; adotv2 = 0. do j=lat0,lat1 do i=lonbeg,lonend im2 = i-2 sinalat = sin(alatm(im2,j)) ! sin(lam) clm2(i,j) = 1.-sinalat*sinalat ! cos^2(lam) be3(i,j) = 1.e-9*be3arr(im2,j) ! be3 is in T (be3arr in nT) do k=lev0,lev1-1 ! ! d_1 = (R_0/R)^1.5 r0or = r0/(r0 + 0.5*(z(k,i,j)+z(k+1,i,j))-h0) rat = 1.e-2*r0or**1.5 ! 1/100 conversion in cm ! ! A_1 dot V = fac( d_1(1) u + d_1(2) v + d_1(3) w adotv1(k,i,j) = rat*( | dvec(im2,j,1,1)*un(k,i,j)+ | dvec(im2,j,2,1)*vn(k,i,j)+ | dvec(im2,j,3,1)*wn(k,i,j)) ! ! Note: clm2 is being used here to represent the squared cosine of the ! quasi-dipole latitude, not of the M(90) latitude, since the wind ! values are aligned vertically, not along the field line. ! rat = rat*sqrt((4.-3.*clm2(i,j))/(4.-3.*r0or*clm2(i,j))) ! ! A_2 dot V = fac( d_2(1) u + d_2(2) v + d_2(3) w adotv2(k,i,j) = rat*( | dvec(im2,j,1,2)*un(k,i,j)+ | dvec(im2,j,2,2)*vn(k,i,j)+ | dvec(im2,j,3,2)*wn(k,i,j)) enddo ! k=lev0,lev1-1 ! ! Calculation of adota(n) = d(n)**2/D ! a1dta2 = (d(1) dot d(2)) /D ! adota1(i,j) = (dvec(im2,j,1,1)**2 + dvec(im2,j,2,1)**2 + | dvec(im2,j,3,1)**2)/dddarr(im2,j) adota2(i,j) = (dvec(im2,j,1,2)**2 + dvec(im2,j,2,2)**2 + | dvec(im2,j,3,2)**2)/dddarr(im2,j) a1dta2(i,j) = (dvec(im2,j,1,1)*dvec(im2,j,1,2) + | dvec(im2,j,2,1)*dvec(im2,j,2,2) + | dvec(im2,j,3,1)*dvec(im2,j,3,2))/dddarr(im2,j) enddo ! i=lon0,lon1 ! ! Must specify levels 1:nlev for addfld because lev0,lev1 are ! actually mlev0,mlev1. ! ! call addfld('Z_ADOTV',' ',' ',z(1:nlev,lon0:lon1,j), ! | 'lev',1,nlev,'lon',lon0,lon1,j) ! ! 3d output fields: ! call addfld('ADOTV1',' ',' ',adotv1(1:nlev,:,j), ! | 'lev',1,nlev,'lon',lon0,lon1,j) ! call addfld('ADOTV2',' ',' ',adotv2(1:nlev,:,j), ! | 'lev',1,nlev,'lon',lon0,lon1,j) enddo ! j=lat0,lat1 ! ! 2d output fields: ! call addfld('ADOTA1',' ',' ',adota1,'lon',lon0,lon1, ! | 'lat',lat0,lat1,0) ! call addfld('ADOTA2',' ',' ',adota2,'lon',lon0,lon1, ! | 'lat',lat0,lat1,0) ! call addfld('A1DTA2',' ',' ',a1dta2,'lon',lon0,lon1, ! | 'lat',lat0,lat1,0) ! call addfld('BE3',' ',' ',be3,'lon',lon0,lon1, ! | 'lat',lat0,lat1,0) end subroutine calc_adotv !----------------------------------------------------------------------- subroutine fieldline_integrals use cons_module,only: ylatm,r0,h0,re ! ! Local: integer :: ier,i,j,k real :: sinlm, clm2, absinim, ra, sqomrra, sqrra, afac, htfac real :: rora,del,omdel,sig1,sig2,ue1,ue2,jgr,jmp,fac,je1pg,je2pg real,dimension(mlon0:mlon1) :: sindm, cosdm, ram, aam, cosiam, | csthdam, rtadram real,dimension(mlon0:mlon1,mlev0:mlev1) :: rrm, sinidm, cosidm, | costhdm, rtramrm, htfunc, htfunc2, je1pg_diag,je2pg_diag ! ! Initialize coefficients: zigm11 = 0. zigm22 = 0. zigm2 = 0. zigmc = 0. rim1 = 0. rim2 = 0. #if defined(INTERCOMM) || defined(CISMAH) zigm1 = 0. #endif je1pg_diag = 0. je2pg_diag = 0. ! ! Subdomain latitude scan for field line integrals: do j=mlat0,mlat1 ! ! Skip only poles (poles will be left 0): ! (If poles are not skipped, clm2 will be zero and it will ! crash at ra = r0/clm2 below) ! if (j==1.or.j==nmlat) cycle ! ! Skip poles and equator: if (j==1.or.j==nmlat.or.j==nmlat/2+1) cycle ! sinlm = sin(ylatm(j)) ! sin(lam_m) clm2 = 1. - sinlm*sinlm ! cos^2(lam_m) absinim = abs(sinlm)/sqrt(1.-.75*clm2) ! | sin I_m | ra = r0/clm2 ! (R_E + H_0)/cos^2(lam_m) sqomrra = sqrt(1.-r0/ra) ! sqrt(1/ R_0/R_A) = sin(lam_m) sqrra = sqrt(r0/ra) ! sqrt(R_0/R_A) afac = 2.*sqrt(ra-r0) ! 2*sqrt(R_A-R_0) (afac is NaN at mag poles) htfac = sqrt(ra-.75*r0) ! sqrt(R_A -3/4*R_0) ! do i=mlon0,mlon1 aam(i) = afac/abs(sini_mag(i,j)) enddo ! ! 2*sqrt( h_A - h_0 )/ |sin I_m | w.r to reference height A(h_R) do k=mlev0,mlev1 do i=mlon0,mlon1 ! ! rr = r0+z-z0 radius of magnetic point ! (Note zpot_mag max was set to h0 in dynamo_inputs) ! rrm(i,k) = r0+zpot_mag(i,j,k)-h0 ! timegcm rrm(i,k) = re+zpot_mag(i,j,k) ! tiegcm ! ! rtramr = ra-r if +ive, zero otherwise rtramrm(i,k) = max(0.,ra-rrm(i,k)) rtramrm(i,k) = sqrt(rtramrm(i,k)) enddo ! i=mlon0,mlon1 enddo ! k=mlev0,mlev1 ! ! Interpolate to midpoints: ! htfunc = factor by which to multiply AAM(I)*d(sqrt(ra-r)) = ds do k=mlev0,mlev1-1 do i=mlon0,mlon1 rrm(i,k) = .5*(rrm(i,k)+rrm(i,k+1)) rtramrm(i,k) = rtramrm(i,k)-rtramrm(i,k+1) htfunc(i,k) = sqrt(ra-.75*rrm(i,k))/htfac htfunc2(i,k) = htfunc(i,k)**2 enddo ! i=mlon0,mlon1 ! write(6,"('fieldline_integrals: j=',i3,' k=',i3,' mlat0,1=', ! | 2i3,' mlon0,1=',2i3,' rrm(:,k)=',/,(6e12.4))") j,k,mlat0, ! | mlat1,mlon0,mlon1,rrm(:,k) ! write(6,"('fieldline_integrals: j=',i3,' k=',i3,' mlat0,1=', ! | 2i3,' mlon0,1=',2i3,' rtramrm(:,k)=',/,(6e12.4))") j,k,mlat0 ! | ,mlat1,mlon0,mlon1,rtramrm(:,k) ! write(6,"('fieldline_integrals: j=',i3,' k=',i3,' mlat0,1=', ! | 2i3,' mlon0,1=',2i3,' htfunc(:,k)=',/,(6e12.4))") j,k,mlat0, ! | mlat1,mlon0,mlon1,htfunc(:,k) ! write(6,"('fieldline_integrals: j=',i3,' k=',i3,' mlat0,1=', ! | 2i3,' mlon0,1=',2i3,' htfunc2(:,k)=',/,(6e12.4))") j,k,mlat0 ! | ,mlat1,mlon0,mlon1,htfunc2(:,k) enddo ! k=mlev0,mlev1 ! ! Compute integrals: do k=mlev0,mlev1-1 do i=mlon0,mlon1 ! ! (R_E+h)/(R_E+h_A) < 1 -> h_A > h rora = min(1.,rrm(i,k)/ra) ! ! (lam_m - lam) / lam_m = ! sqrt(1-r_0/r_A)sqrt(r/r_A) - sqrt(r_0/r_A)sqrt(1-r/r_A) del = (sqomrra*sqrt(rora)-sqrra*sqrt(1.-rora))/ | abs(ylatm(j)) omdel = 1. - del ! ! Interpolate conductivities and winds in latitude along field line, assuming ! linear variation between foot of field line and magnetic equator. ! (For field lines other than those near the magnetic equator, del is nearly ! zero, so that the interpolated values are essentially the values for the ! latitude of the foot of the field line; inaccuracy of the assumption of ! linear variation is thus unimportant for these field lines.) ! ! Here, mag equator ped_meq, etc. are from mp_mageq, called from dynamo inputs: sig1 = omdel*ped_mag(i,j,k) + del*ped_meq(i,k) sig2 = omdel*hal_mag(i,j,k) + del*hal_meq(i,k) ue1 = omdel*adotv1_mag(i,j,k) + del*adotv1_meq(i,k) ue2 = omdel*adotv2_mag(i,j,k) + del*adotv2_meq(i,k) ! write(6,"('fieldline_integrals: i=',i3,' k=',i3,' j=',i3, ! | ' omdel=',e12.4,' del=',e12.4,' sig1,2=',2e12.4,' ue1,2=', ! | 2e12.4)") i,k,j,omdel,del,sig1,sig2,ue1,ue2 ! ! height varying factors: ds = aam*htfunc ! d_1^2/D = 1/htfunc * adota1_mag(i,j) ! d_2^2/D = htfunc * adota2_mag(i,j) ! d_1*d_2/D = 1 * a1dta2_mag(i,j) ! ! zigm11: int (sigma_p*d_1^2/D) ds : d_1^2/D ! zigm22: int (sigma_p*d_2^2/D) ds : d_2^2/D ! zigm11(i,j) = zigm11(i,j) + sig1*rtramrm(i,k) zigm22(i,j) = zigm22(i,j) + sig1*rtramrm(i,k)*htfunc2(i,k) ! ! zigmc: int (sigma_p*d_1*d_2/D) ds ! zigm2: int (sigma_h) ds ! zigmc(i,j) = zigmc(i,j) + sig1*rtramrm(i,k)*htfunc(i,k) zigm2(i,j) = zigm2(i,j) + sig2*rtramrm(i,k)*htfunc(i,k) #if defined(INTERCOMM) || defined(CISMAH) zigm1(i,j) = zigm1(i,j) + sig1*rtramrm(i,k)*htfunc(i,k) #endif ! ! rim1: int [sigma_p*d_1^2/D u_e2+(sigma_h-(sigma_p*d_1*d_2)/D) u_e1] ds ! rim2: int [(sigma_h+sigma_p*d_1*d_2/D) u_e2-sigma_p*d_2^2/D u_e1 ] ds ! rim1(i,j) = rim1(i,j) + (sig1*adota1_mag(i,j)*ue2 + | (sig2 - sig1*a1dta2_mag(i,j))*htfunc(i,k)*ue1)*rtramrm(i,k) rim2(i,j) = rim2(i,j) + (sig1*adota2_mag(i,j)*htfunc2(i,k)* | ue1 - (sig2 + sig1*a1dta2_mag(i,j))*htfunc(i,k)*ue2)* | rtramrm(i,k) ! ! rim1 add = 1/B_e3 * int [-k(ti+te)*grad Ne + rho*g]xB / B^2*d1/D ds ! rim2 add = 1/B_e3 * int [-k(ti+te)*grad Ne + rho*g]xB / B^2*d2/D ds ! je1pg & je2pg [A/cm^2] ! if (current_pg > 0) then fac = be3_mag(i,j) fac = rtramrm(i,k)*htfunc(i,k)/fac ! ds/be3/aam je1pg = omdel*je1pg_mag(i,j,k) + del*je1pg_meq(i,k) je2pg = omdel*je2pg_mag(i,j,k) + del*je2pg_meq(i,k) rim1(i,j) = rim1(i,j)+ 1.e4*je1pg*fac ! 1e4 for A/cm^2 -> A/m^2 rim2(i,j) = rim2(i,j)+ 1.e4*je2pg*fac ! 1e4 for A/cm^2 -> A/m^2 endif enddo ! i=mlon0,mlon1 enddo ! k=mlev0,mlev1 ! ! call addfld('ZIGM11_a','ZIGM11_a integrals before *.01','S', ! | zigm11(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('ZIGM22_a','ZIGM22_a integrals before *.01','S', ! | zigm22(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('ZIGMC_a','ZIGMC_a integrals before *.01','S', ! | zigmc(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('ZIGM2_a','ZIGM2_a integrals before *.01','S', ! | zigm2(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('RIM1_a','RIM1_a integrals before *.01','A/m', ! | rim1(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('RIM2_a','RIM2_a integrals before *.01','A/m', ! | rim2(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! ! Complete calculation and place result in /coefm/ zigm's ! rim's are in A/m multiply by 1/100 to convert from [cm] to [m] ! ! At this point, ! zigm11 is int[sig_p*d_1^2/D] ds, i.e. Sigma_(phi phi)/abs(sin Im) ! zigm22 is int[sig_p*d_2^2/D] ds, i.e. Sigma_(lam lam)*abs(sin Im) ! zigmc is int[sig_p*d_1*d_2/D] ds, i.e. Sigma_c ! zigm2 is int[sigma_h] ds, i.e. Sigma_h ! ! rim1 is int[(sigma_h-sigma_p*d_1*d_2/D)u_e1 + sigma_p*d_1^2/D u_e2] *A(h_r)* ! B_e3 ds, i.e. K_(m phi)^D/abs(sin Im) ! rim2 is int[(sigma_h+sigma_p*d_1*d_2/D)u_e2 - sigma_p*d_2^2/D u_e1] *A(h_r)* ! B_e3 ds, K_(m lam)^D ( minus in northern hemisphere ! Change sign of RIM(2) in S. hemisphere to be compatible with transf ! At this point, rim2 is +-K_(m lam)^D ! do i = mlon0,mlon1 zigm11(i,j) = 1.e-2*zigm11(i,j)*aam(i)*adota1_mag(i,j) zigm22(i,j) = 1.e-2*zigm22(i,j)*aam(i)*adota2_mag(i,j) zigmc(i,j) = 1.e-2*zigmc (i,j)*aam(i)*a1dta2_mag(i,j) zigm2(i,j) = 1.e-2*zigm2 (i,j)*aam(i) #if defined(INTERCOMM) || defined(CISMAH) zigm1(i,j) = 1.e-2*zigm1(i,j)*aam(i) #endif rim1(i,j) = 1.e-2*rim1(i,j)*aam(i)*be3_mag(i,j) rim2(i,j) = 1.e-2*rim2(i,j)*aam(i)*be3_mag(i,j) enddo ! i = 1,nmlon ! ! 12/9/11 btf: Similiar diffs w/ odyn as above, with addition of ! some diffs apparently in geo pole positions (presumably from ! diffs in adota1,2,be3) Overall structure is ok, diffs not visible ! in full-field plots. ! ! call addfld('ZIGM11_b','ZIGM11_b integrals after *.01','S', ! | zigm11(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('ZIGM22_b','ZIGM22_b integrals after *.01','S', ! | zigm22(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('ZIGMC_b','ZIGMC_b integrals after *.01','S', ! | zigmc(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('ZIGM2_b','ZIGM2_b integrals after *.01','S', ! | zigm2(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('RIM1_b','RIM1_b integrals after *.01','A/m', ! | rim1(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) ! call addfld('RIM2_b','RIM2_b integrals after *.01','A/m', ! | rim2(mlon0:mlon1,j),'mlon',mlon0,mlon1,'mlat',j,j,0) enddo ! j=mlat0,mlat1 (without poles) ! end subroutine fieldline_integrals !----------------------------------------------------------------------- subroutine alloc_pdyn ! ! Allocate and initialize arrays for parallel dynamo (module data) ! (called once per run from allocdata.F) ! integer :: istat,n integer :: mlon00,mlon11,mlat00,mlat11 ! mlon00=mlon0-1 ; mlon11=mlon1+1 mlat00=mlat0-1 ; mlat11=mlat1+1 ! ! 2d fields on mag subdomains (i,j): ! Certain fields are allocated with halos mlon0-1:mlon1+1,mlat0-1:mlat1+1 ! allocate(zigm11(mlon00:mlon11,mlat00:mlat11),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: zigm11') zigm11 = 0. allocate(zigmc(mlon00:mlon11,mlat00:mlat11) ,stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: zigmc') zigmc = 0. allocate(zigm2(mlon00:mlon11,mlat00:mlat11) ,stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: zigm2') zigm2 = 0. allocate(zigm22(mlon00:mlon11,mlat00:mlat11),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: zigm22') zigm22 = 0. allocate(rhs(mlon00:mlon11,mlat00:mlat11) ,stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: rhs') rhs = 0. allocate(rim1(mlon00:mlon11,mlat00:mlat11) ,stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: rim1') rim1 = 0. allocate(rim2(mlon00:mlon11,mlat00:mlat11) ,stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: rim2') rim2 = 0. allocate(phimsolv(mlon00:mlon11,mlat00:mlat11),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: phimsolv') phimsolv = 0. allocate(phim2d(mlon00:mlon11,mlat00:mlat11),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: phim2d') phim2d = 0. #if defined(INTERCOMM) || defined(CISMAH) allocate(zigm1(mlon00:mlon11,mlat00:mlat11),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: zigm1') zigm1 = 0. allocate(nsrhs(mlon00:mlon11,mlat00:mlat11),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: nsrhs') nsrhs = 0. #endif ! ! 3d phim and electric field on mag subdomains: allocate(phim3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: phim3d') phim3d = 0. allocate(ed13d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: ed13d') ed13d = 0. allocate(ed23d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: ed23d') ed23d = 0. allocate(ephi3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: ephi3d') ephi3d = 0. allocate(elam3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: elam3d') elam3d = 0. allocate(emz3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: emz3d') emz3d = 0. allocate(zpotm3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: zpotm3d') zpotm3d = 0. ! ! 3d electric field on geographic subdomains: allocate(ex(nlevp1,lon0-2:lon1+2,lat0:lat1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: ex') allocate(ey(nlevp1,lon0-2:lon1+2,lat0:lat1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: ey') allocate(ez(nlevp1,lon0-2:lon1+2,lat0:lat1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: ez') ex=0. ; ey=0. ; ez=0. ! ! 3d electric potential on geographic subdomains (regridded from phim3d): allocate(phig3d(lon0:lon1,lat0:lat1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: phig3d') phig3d = 0. allocate(coef_rhs(mlon0:mlon1,mlat0:mlat1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: coef_rhs') coef_rhs = 0. allocate(coef_cofum(mlon0:mlon1,mlat0:mlat1,1:9),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: coef_cofum') coef_cofum = 0. ! allocate(ped_meq(mlon0:mlon1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: ped_meq') ped_meq = 0. allocate(hal_meq(mlon0:mlon1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: hal_meq') hal_meq = 0. allocate(adotv1_meq(mlon0:mlon1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: adotv1_meq') adotv1_meq = 0. allocate(adotv2_meq(mlon0:mlon1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: adotv2_meq') adotv2_meq = 0. allocate(je1pg_meq(mlon0:mlon1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: je1pg_meq') je1pg_meq = 0. allocate(je2pg_meq(mlon0:mlon1,mlev0:mlev1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: je2pg_meq') je2pg_meq = 0. ! ! Fields input to mp_mageq (6 fields at full mag subdomain i,j,k): write(6,"('alloc_pdyn: allocate fmeq_in with mlev0,1=',2i4)") | mlev0,mlev1 allocate(fmeq_in(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1,6), | stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: fmeq_in') fmeq_in = 0. ! ! Fields output by mp_mageq (6 fields at mag subdomain i,k) allocate(fmeq_out(mlon0:mlon1,mlev0:mlev1,6),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: fmeq_out') fmeq_out = 0. ! am 9/14/2010 ! for L2norm (only for checking not necessary and can be removed later) allocate(ressolv(mlon0:mlon1,mlat0:mlat1),stat=istat) if (istat /= 0) call shutdown('alloc_pdyn: ressolv') ressolv = 0. end subroutine alloc_pdyn !----------------------------------------------------------------------- subroutine complete_integrals ! ! Field line integrals for each hemisphere have been calculated in ! mag subdomains. Now, complete these arrays with equator and polar ! values, and sum the integrals from the 2 hemispheres. ! (most of this was in sub transf in old serial dynamo) ! This is done by obtaining global mag fields via mpi exchange ! of mag subdomains, completing the global fields, and updating ! the subdomains. ! The 6 2d fields are: zigm11,zigm22,zigmc,zigm2,rim1,rim2 ! use cons_module,only: rcos0s,dt1dts use mpi_module,only: mp_mageq_jpm1,mp_magpole_2d,mp_mag_foldhem use mpi_module,only: mytid #if defined(INTERCOMM) || defined(CISMAH) use mpi_module,only: mp_mag_halos use cism_coupling_module,only: gzigm1,gzigm2,gnsrhs,cism_ucurrent use my_esmf,only: | mag_zigm1, mag_zigm2, mag_nsrhs, | geo_zigm1, geo_zigm2, geo_nsrhs #endif implicit none ! ! Local: integer :: i,ii,j,ifld,lonend real :: fmsub(mlon0:mlon1,mlat0:mlat1,nf2d) real :: corfac real,parameter :: unitvm(nmlon)=1. ! ! External: real,external :: sddot ! util.F ! ! For equatorial values, we need latitudes eq+1 and eq-1: ! Local feq_jpm1(nmlonp1,2,6) is returned by mp_mageq_jpm1, ! where the 2 dim contains lats nmlath-1, nmlath+1. These ! are global in lon, even tho each subd uses only its own i's. ! These mag equator values do not show up on plots because ! of the small factor .06 and .125. ! The 6 2d fields are: zigm11,zigm22,zigmc,zigm2,rim1,rim2 ! Must specify mlon0:mlon1,mlat0:mlat1 because these fields ! were allocated to include single-point halos (these calls ! exclude the halo points): ! fmsub(:,:,1) = zigm11(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,2) = zigm22(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,3) = zigmc (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,4) = zigm2 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,5) = rim1 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1) #if defined(INTERCOMM) || defined(CISMAH) fmsub(:,:,7) = zigm1(mlon0:mlon1,mlat0:mlat1) #endif call mp_mageq_jpm1(fmsub,mlon0,mlon1,mlat0,mlat1,nmlonp1, | feq_jpm1,nf2d) ! ! From sub fieldline_integrals: ! zigm11 is int[sig_p*d_1^2/D] ds, i.e. Sigma_(phi phi)/abs(sin Im) ! zigm22 is int[sig_p*d_2^2/D] ds, i.e. Sigma_(lam lam)*abs(sin Im) ! zigmc is int[sig_p*d_1*d_2/D] ds, i.e. Sigma_c ! zigm2 is int[sigma_h] ds, i.e. Sigma_h ! ! rim1 is int[(sigma_h-sigma_p*d_1*d_2/D)u_e1 + sigma_p*d_1^2/D u_e2] *A(h_r)* ! B_e3 ds, i.e. K_(m phi)^D/abs(sin Im) ! rim2 is int[(sigma_h+sigma_p*d_1*d_2/D)u_e2 - sigma_p*d_2^2/D u_e1] *A(h_r)* ! B_e3 ds, K_(m lam)^D ( minus in northern hemisphere ! Change sign of RIM(2) in S. hemisphere to be compatible with transf ! At this point, rim2 is +-K_(m lam)^D ! ! Equatorial values: ! Assume that quantities primarily dependent on Pedersen conductivity ! have field-line integrals 1/4 as large as the averages for next-higher ! field lines; quantities primarily dependent on Hall conductivity ! have field-line integrals 0.12 as large as the averages for next-higher ! field lines. Exact values chosen should not be important for potential ! calculation, as long as they are physically reasonable and not too ! different from adjacent values. ! do j=mlat0,mlat1 if (j==jlateq) then zigm11(mlon0:mlon1,j) = .125*(feq_jpm1(mlon0:mlon1,1,1)+ | feq_jpm1(mlon0:mlon1,2,1)) zigm22(mlon0:mlon1,j) = .125*(feq_jpm1(mlon0:mlon1,1,2)+ | feq_jpm1(mlon0:mlon1,2,2)) zigmc (mlon0:mlon1,j) = .125*(feq_jpm1(mlon0:mlon1,1,3)+ | feq_jpm1(mlon0:mlon1,2,3)) zigm2 (mlon0:mlon1,j) = .060*(feq_jpm1(mlon0:mlon1,1,4)+ | feq_jpm1(mlon0:mlon1,2,4)) rim1 (mlon0:mlon1,j) = .060*(feq_jpm1(mlon0:mlon1,1,5)+ | feq_jpm1(mlon0:mlon1,2,5)) rim2 (mlon0:mlon1,j) = .060*(feq_jpm1(mlon0:mlon1,1,6)+ | feq_jpm1(mlon0:mlon1,2,6)) #if defined(INTERCOMM) || defined(CISMAH) zigm1(mlon0:mlon1,j) = .060*(feq_jpm1(mlon0:mlon1,1,7)+ | feq_jpm1(mlon0:mlon1,2,7)) #endif ! ! Include the boundary condition at the equator eq.(5.30) in ! Richmond (1995) Ionospheric Electrodynamics use. Mag. Apex Coord. ! J.Goemag.Geoelectr. 47,191-212 ! Sig_phiphi/abs(sin Im) = 0.5*Sig_cowling/abs(sin Im) ! = 0.5/abs(sin Im)*(Sig_phiphi - Sig_philam*sig_lamphi/Sig_lamlam) ! = 0.5/abs(sin Im)*(Sig_phiphi + (Sig_h-sig_c)*(Sig_h+sig_c)/Sig_lamlam) ! rim(1) / |sin I_m| = I_1 = R/2*(K_mphi - Sig_philam/Sig_lamlam*K_mlam) ! do i=mlon0,mlon1 zigm11(i,j) = zigm11(i,j)+(zigm2(i,j)-zigmc(i,j))* | (zigm2(i,j)+zigmc(i,j))/zigm22(i,j) rim1(i,j) = rim1(i,j) - (zigm2(i,j)-zigmc(i,j))/ | zigm22(i,j)*rim2(i,j) ! ! 12/9/11 btf: current dynamo.F (tiegcm and timegcm) has the following ! 2 "do-nothing" lines, but earlier versions (e.g., tag timegcm1-2dev7 ! has 0.5* for both lines, like the commented lines below. ! zigm11(i,j) = zigm11(i,j) ! current dynamo.F rim1(i,j) = rim1(i,j) ! current dynamo.F ! zigm11(i,j) = 0.5*zigm11(i,j) ! old timegcm1-2dev7 tag ! rim1(i,j) = 0.5*rim1(i,j) ! old timegcm1-2dev7 tag enddo ! i=mlon0,mlon1 endif ! j at equator enddo ! j=mlat0,mlat1 ! ! zigm11, etc are allocated to include halo points, so must ! specify subdomains subsections here: ! ! call addfld('ZIGM11_c','ZIGM11_c after mageq','S', ! | zigm11(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGM22_c','ZIGM22_c after mageq','S', ! | zigm22(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGMC_c','ZIGMC_c after mageq','S', ! | zigmc(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGM2_c','ZIGM2_c after mageq','S', ! | zigm2(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('RIM1_c','RIM1_c after mageq','A/m', ! | rim1(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('RIM2_c','RIM2_c after mageq','A/m', ! | rim2(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! ! Using notation of Richmond (1995) on right-hand side below: ! Sigma_(phi phi) = zigm11*abs(sin I_m) ! Sigma_(lam lam) = zigm22/abs(sin I_m) ! Sigma_(phi lam) = +-(zigm2-zigmc) ! Sigma_(lam phi) = -+(zigm2+zigmc) ! K_(m phi)^D = rim(1)*abs(sin I_m) ! K_(m lam)^D = +-rim(2) ! ! Transforming PDE from original apex (theta_a) to new apex grid (theta_0) ! which is equally spaced in magnetic latitude ! SCALE quantities to modified (0) magnetic latitude system, multiplying or ! dividing by abs(sin I_m) [inverse contained in DT1DTS] as necessary. ! Sign of K_(m lam)^D in southern hemisphere remains reversed. ! for the mixed terms the transformation from the integration and differentiation ! canceled out (zigmc, zigm2) ! DT1DTS : d theta_0/ d theta_a / abs(sin I_m) ! RCOS0S : cos(theta_0)/ cos(theta_a) ! ! corfac: abs(I_m)*d theta_a/d theta_0 * cos(theta_0)/ cos(theta_a) ! zigm11: abs(I_m)*d theta_a/d theta_0 * cos(theta_0)/ cos(theta_a) ! zigm22: 1/abs(I_m)*d theta_0/d theta_a * cos(theta_a)/ cos(theta_0) ! rim(1): abs(I_m)*d theta_a/d theta_0 ! rim(2): cos(theta_a)/ cos(theta_0) ! do j=mlat0,mlat1 if (j==1.or.j==nmlat) cycle ! skip poles corfac = rcos0s(j)/dt1dts(j) do i=mlon0,mlon1 zigm11(i,j) = zigm11(i,j)*corfac zigm22(i,j) = zigm22(i,j)/corfac rim1(i,j) = rim1(i,j)/dt1dts(j) rim2(i,j) = rim2(i,j)/rcos0s(j) enddo enddo ! ! For polar values, we need south pole plus 1 and 2 (j==2,3), ! and north pole minus 1 and 2 (j==nmlat-1,nmlat-2). These ! are returned by sub mp_magpole_jpm2 (mpi.F): ! Must specify (mlon0:mlon1,mlat0:mlat1) because zigmxx and rims ! are allocated to include halo cells. ! fmsub(:,:,1) = zigm11(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,2) = zigm22(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,3) = zigmc (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,4) = zigm2 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,5) = rim1 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1) #if defined(INTERCOMM) || defined(CISMAH) fmsub(:,:,7) = zigm1 (mlon0:mlon1,mlat0:mlat1) #endif ! ! mp_magpole_2d returns fpole_jpm2(nmlonp1,1->4,nf) as: ! 1: j = 2 (spole+1) ! 2: j = 3 (spole+2) ! 3: j = nmlat-1 (npole-1) ! 4: j = nmlat-2 (npole-2) ! ! subroutine mp_magpole_2d(f,ilon0,ilon1,ilat0,ilat1, ! nglblon,jspole,jnpole,fpole_jpm2,nf) call mp_magpole_2d(fmsub,mlon0,mlon1,mlat0,mlat1,nmlonp1, | 1,nmlat,fpole_jpm2,nf2d) ! write(6,"('complete_integrals: zigm11 global lons at j==2 = ',/, ! | (6e12.4))") fpole_jpm2(:,1,1) ! write(6,"('complete_integrals: zigm11 global lons at j==nmlat-1 =' ! | ,/,(6e12.4))") fpole_jpm2(:,3,1) ! write(6,"('complete_integrals: zigm22 global lons at j==2 = ',/, ! | (6e12.4))") fpole_jpm2(:,1,2) ! write(6,"('complete_integrals: zigm22 global lons at j==nmlat-1 =' ! | ,/,(6e12.4))") fpole_jpm2(:,3,2) ! ! the PDE is divided by 1/ DT0DTS ! Sigma_(phi phi) = zigm11/ rcos0s * dt0dts ! Sigma_(lam lam) = zigm22 * rcos0s / dt0dts ! Sigma_(phi lam) = +-(zigm2-zigmc) ! Sigma_(lam phi) = -+(zigm2+zigmc) ! K_(m phi)^D = rim(1) * dt0dts ! K_(m lam)^D = +-rim(2) * rcos0s ! ! Compute polar values for the conductances, 4th order interpolation: ! ! do j=mlat0,mlat1 ! ! South pole: if (j==1) then ! south pole (use fpole_jpm2(nmlon,1->2,nf) zigm11(mlon0,j)=(4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,1,1))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,2,1)))/ | (3.*float(nmlon)) zigm22(mlon0,j)=(4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,1,2))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,2,2)))/ | (3.*float(nmlon)) zigmc (mlon0,j)=(4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,1,3))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,2,3)))/ | (3.*float(nmlon)) zigm2 (mlon0,j)=(4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,1,4))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,2,4)))/ | (3.*float(nmlon)) #if defined(INTERCOMM) || defined(CISMAH) zigm1(mlon0,j) = (4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,1,7))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,2,7)))/ | (3.*float(nmlon)) #endif ! ! Extend south pole over longitude: do i=mlon0+1,mlon1 zigm11(i,j) = zigm11(mlon0,j) zigm22(i,j) = zigm22(mlon0,j) zigmc (i,j) = zigmc (mlon0,j) zigm2 (i,j) = zigm2 (mlon0,j) #if defined(INTERCOMM) || defined(CISMAH) zigm1 (i,j) = zigm1 (mlon0,j) #endif enddo ! i=mlon0,mlon1 ! ! RHS vector (I_1,I_2): average over south pole: ! (use fpole_jpm2(i,1,nf), i.e. j==2, and lons across the pole) lonend = mlon1 if (mlon1==nmlonp1) lonend = mlon1-1 do i=mlon0,lonend ii = 1+mod(i-1+nmlon/2,nmlon) rim1(i,j) = 0.5*(fpole_jpm2(i,1,5)-fpole_jpm2(ii,1,5)) rim2(i,j) = 0.5*(fpole_jpm2(i,1,6)-fpole_jpm2(ii,1,6)) enddo ! ! North pole: elseif (j==nmlat) then ! north pole (use fpole_jpm2(nmlon,3->4,1,nf) zigm11(mlon0,j)=(4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,3,1))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,4,1)))/ | (3.*float(nmlon)) zigm22(mlon0,j)=(4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,3,2))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,4,2)))/ | (3.*float(nmlon)) zigmc (mlon0,j)=(4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,3,3))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,4,3)))/ | (3.*float(nmlon)) zigm2 (mlon0,j)=(4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,3,4))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,4,4)))/ | (3.*float(nmlon)) #if defined(INTERCOMM) || defined(CISMAH) zigm1(mlon0,j) = (4.* | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,3,7))- | sddot(nmlon,unitvm,fpole_jpm2(1:nmlon,4,7)))/ | (3.*float(nmlon)) #endif ! ! Extend north pole over longitude: do i=mlon0+1,mlon1 zigm11(i,j) = zigm11(mlon0,j) zigm22(i,j) = zigm22(mlon0,j) zigmc (i,j) = zigmc (mlon0,j) zigm2 (i,j) = zigm2 (mlon0,j) #if defined(INTERCOMM) || defined(CISMAH) zigm1 (i,j) = zigm1 (mlon0,j) #endif enddo ! i=mlon0,mlon1 ! ! RHS vector (I_1,I_2): average over north pole: ! (use fpole_jpm2(i,3,nf), i.e. j==nmlat-1, and lons across the pole) lonend = mlon1 if (mlon1==nmlonp1) lonend = mlon1-1 do i=mlon0,lonend ii = 1+mod(i-1+nmlon/2,nmlon) rim1(i,j) = 0.5*(fpole_jpm2(i,3,5)-fpole_jpm2(ii,3,5)) rim2(i,j) = 0.5*(fpole_jpm2(i,3,6)-fpole_jpm2(ii,3,6)) enddo endif ! south or north pole enddo ! j=mlat0,mlat1 ! ! call addfld('ZIGM11_d','ZIGM11_d after mag poles','S', ! | zigm11(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGM22_d','ZIGM22_d after mag poles','S', ! | zigm22(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGMC_d','ZIGMC_d after mag poles','S', ! | zigmc(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGM2_d','ZIGM2_d after mag poles','S', ! | zigm2(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('RIM1_d','RIM1_d after mag poles','A/m', ! | rim1(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('RIM2_d','RIM2_d after mag poles','A/m', ! | rim2(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) #if defined(INTERCOMM) || defined(CISMAH) ! ! ! Calculate global neutral winds generated current to transfer to CISM M-I couplier ! call mp_mag_halos(rim1,mlon0,mlon1,mlat0,mlat1,1,1) call mp_mag_halos(rim2,mlon0,mlon1,mlat0,mlat1,1,1) call cism_ucurrent(rim1,rim2,mlon0,mlon1,mlat0,mlat1,nsrhs) ! call cism_ucurrent ! ! Transform height-integrated parameters from APEX coordinates to geographic ! coordinates to be passed to the M-I coupled in the CISM ! call mag2geo_2d(zigm1(mlon0:mlon1,mlat0:mlat1),gzigm1,mag_zigm1, | geo_zigm1,'ZIGM1 ') call mag2geo_2d(zigm2(mlon0:mlon1,mlat0:mlat1),gzigm2,mag_zigm2, | geo_zigm2,'ZIGM2 ') call mag2geo_2d(nsrhs(mlon0:mlon1,mlat0:mlat1),gnsrhs,mag_nsrhs, | geo_nsrhs,'NSRHS ') #endif ! ! Set coefficients for current density calculations. ! Note nosocoef is a stub in current.F90 (but not in the ! current module) which calls noso_coef. This is necessary ! to avoid circular dependency between the pdynamo and current ! modules. Subs noso_crrt and noso_crdens (current.F90) are ! called from advance after pdynamo. ! if (current_kq > 0) then call gather_pdyn if (mytid==0) call nosocoef endif ! ! Fold south hemisphere over onto north, and set periodic points: fmsub(:,:,1) = zigm11(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,2) = zigm22(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,3) = zigmc (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,4) = zigm2 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,5) = rim1 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1) #if defined(INTERCOMM) || defined(CISMAH) fmsub(:,:,7) = zigm1 (mlon0:mlon1,mlat0:mlat1) #endif call mp_mag_foldhem(fmsub,mlon0,mlon1,mlat0,mlat1,nf2d) call mp_mag_periodic_f2d(fmsub,mlon0,mlon1,mlat0,mlat1,nf2d) zigm11(mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,1) zigm22(mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,2) zigmc (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,3) zigm2 (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,4) rim1 (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,5) rim2 (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,6) #if defined(INTERCOMM) || defined(CISMAH) zigm1 (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,7) #endif ! ! Reverse sign of zigmc in northern hemisphere. do j=mlat0,mlat1 if (j >= nmlath) then zigmc(mlon0:mlon1,j) = -zigmc(mlon0:mlon1,j) endif enddo ! do j=mlat0,mlat1 ! if (j==nmlath) write(6,"('j=nmlath=',i3,' mlon0,1=',2i3, ! | ' zigmc(mlon0:mlon1,j) after sign reversal =',/,(6e12.4))") ! | j,mlon0,mlon1,zigmc(mlon0:mlon1,j) ! enddo ! call addfld('ZIGM11_e','ZIGM11_e after fold hem','S', ! | zigm11(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGM22_e','ZIGM22_e after fold hem','S', ! | zigm22(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGMC_e','ZIGMC_e after fold hem','S', ! | zigmc(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('ZIGM2_e','ZIGM2_e after fold hem','S', ! | zigm2(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('RIM1_e','RIM1_e after fold hem','A/m', ! | rim1(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! call addfld('RIM2_e','RIM2_e after fold hem','A/m', ! | rim2(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) ! end subroutine complete_integrals !----------------------------------------------------------------------- subroutine rhspde use cons_module,only: pi_dyn,dlatm,dlonm,r0 use cons_module,only: rcos0s,dt1dts ! ! Calculate right-hand side from rim1,2 on mag subdomains. ! Use global longitude arrays for poles and equator obtained ! by sub complete_integrals. ! ! Local: integer :: j,i,ii real :: pi real,dimension(nmlat) :: tint1 real,parameter :: unitvm(nmlon)=1. real :: rim2_npm1(nmlonp1), ! global rim2 at nmlat-1 | rim2_eqp1(nmlonp1), ! global rim2 at meq+1 | rim2_meq (nmlonp1), ! global rim2 at mag eq | rim2_tmp (nmlonp1), ! tmp array | rim1_meq (nmlonp1), ! global rim1 at mag eq | zigm2_meq(nmlonp1), ! needed for rim1_meq | zigmc_meq(nmlonp1), ! needed for rim1_meq | zigm22_meq(nmlonp1) ! needed for rim1_meq ! real :: fsave(mlon0:mlon1,mlat0:mlat1) ! ! External: real,external :: sddot ! util.F ! ! call addfld('RIM1',' ','S',rim1(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('RIM2',' ','S',rim2(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) pi = pi_dyn do j=1,nmlat tint1(j) = cos(-pi/2.+(j-1)*dlatm) enddo ! ! Init rhs subdomains: rhs(:,:) = 0. ! ! Will need rim2 at npole-1 and mag equator: ! rim2_npm1: global rim2 at nmlat-1: rim2_npm1(:) = fpole_jpm2(:,3,6)+fpole_jpm2(:,1,6) ! rim2_meq: global rim2 at mag equator jlateq: rim2_meq(:) = .060*(feq_jpm1(:,1,6)+feq_jpm1(:,2,6)) rim2_meq(:) = rim2_meq(:)/rcos0s(jlateq) rim2_meq(:) = rim2_meq(:)*2. ! fold eq on itself ! ! Perform differentiation of rim(2) w.r.t. lam_0: ! +/- (d [ K_(m lam)^D * cos(lam_m)]/ d lam_0 ) /cos ( lam_0) = ! + (d [ K_(m lam)^D * cos(lam_m)]/ d lam_m ) /cos ( lam_m) / (RCOS0S*DT0DTS) = ! +/- (d [ K_(m lam)^D(0) * cos(lam_0)]/ d lam_0 ) /cos ( lam_0) = ! ! Lat scan to define rhs subdomains: do j=mlat0,mlat1 ! ! North Pole (redundant in longitude): if (j == nmlat) then ! north pole do i=mlon0,mlon1 rhs(i,j) = -2./float(nmlon)*sddot(nmlon,unitvm,rim2_npm1)/ | tint1(nmlat-1) enddo ! ! Include the boundary condition at the equator. ! rhs(equator)/R = d (K_mphi^DT(0) - sig_philam/sig_lamlam*K_mlam^DT(0)) / d phi_m ! + d (cos lam_0 * K_mlam^DT(0))/ d lam_0 ! from Cicely's notes: ! I_1 = 0.5*(K_(m phi)^DT(0) - Sig_(phi lam)/Sig_(lam lam)*K_(ml am)^DT(0)) ! I_2 = K_(m lam)^DT(0) ! differentiate ! rhs = (I_1(i+1/2,j)-I_1(i-1/2,j))/dlonm + ! (2*cos(lam_0)_(j+1/2)*I_2(i,j+1/2))/dlat_0 ! (first calc global mag equator as in complete_integrals) ! elseif (j == jlateq) then ! mag equator rim2_eqp1(:) = feq_jpm1(:,2,6)/rcos0s(j-1)+ | feq_jpm1(:,1,6)/rcos0s(j+1) zigm22_meq(:) = .125*(feq_jpm1(:,1,2)+feq_jpm1(:,2,2)) zigmc_meq (:) = .125*(feq_jpm1(:,1,3)+feq_jpm1(:,2,3)) zigm2_meq (:) = .060*(feq_jpm1(:,1,4)+feq_jpm1(:,2,4)) rim1_meq (:) = .060*(feq_jpm1(:,1,5)+feq_jpm1(:,2,5)) rim2_tmp (:) = .060*(feq_jpm1(:,1,6)+feq_jpm1(:,2,6)) do i=1,nmlonp1 rim1_meq(i) = rim1_meq(i) - (zigm2_meq(i)-zigmc_meq(i))/ | zigm22_meq(i)*rim2_tmp(i) enddo rim1_meq(:) = rim1_meq(:)/dt1dts(j) rim1_meq(:) = rim1_meq(:)*2. ! fold eq on itself do i=mlon0,mlon1 if (i==1) then ! western most lon rhs(i,j) = 0.5/dlonm*(rim1(i+1,j)-rim1_meq(nmlon)) rhs(i,j) = rhs(i,j)+1./dlatm*(tint1(j)*rim2(i,j)+ | tint1(j+1)*rim2_eqp1(i)) elseif (i==nmlonp1) then ! eastern most lon rhs(i,j) = 0.5/dlonm*(rim1_meq(1)-rim1(i-1,j)) rhs(i,j) = rhs(i,j)+1./dlatm*(tint1(j)*rim2(i,j)+ | tint1(j+1)*rim2_eqp1(i)) else ! body of i subdomain ! Note that rim1 halos were set before calling this subroutine. rhs(i,j) = 0.5/dlonm*(rim1(i+1,j)-rim1(i-1,j)) rhs(i,j) = rhs(i,j)+1./dlatm*(tint1(j)*rim2(i,j)+ | tint1(j+1)*rim2_eqp1(i)) endif enddo ! i=mlon0,mlon1 ! write(6,"('rhs: j==jlateq=',i4,' rhs(mlon0:mlon1,j)=',/, ! | (6e12.4))") j,rhs(mlon0:mlon1,j) ! ! North hemisphere (not npole and not equator): ! (allow south hemisphere to remain 0) ! (use rim1 instead of tint33) ! elseif (j > jlateq) then ! north hem only (excluding npole) do i=mlon0,mlon1 rhs(i,j) = 1./(dlonm*tint1(j))*0.5*(rim1(i+1,j)-rim1(i-1,j)) if (j == jlateq+1) then rhs(i,j) = rhs(i,j)+1./(dlatm*tint1(j))* | 0.5*(rim2(i,j+1)*tint1(j+1)-rim2_meq(i)*tint1(j-1)) else rhs(i,j) = rhs(i,j)+1./(dlatm*tint1(j))* | 0.5*(rim2(i,j+1)*tint1(j+1)-rim2(i,j-1)*tint1(j-1)) endif enddo ! write(6,"('rhs: j=nhem=',i4,' rhs(mlon0:mlon1,j)=',/, ! | (6e12.4))") j,rhs(mlon0:mlon1,j) endif enddo ! j=mlat0,mlat1 ! ! scale (multiply by earth radius in meter = R0*1.E-2) ![( d K_(m phi)^D / d phi /(cos(theta_m)?) + ! (d [ K_(m lam)^D * cos(lam_m)]/ d lam_m ) /cos ( lam_m) ] * R / (RCOS0S*DT0DTS) ! ~ J_(Mr)*r^2*cos(theta_m)/cos(theta_0)/DT0DTS ! theta_m = theta_s ! do j=mlat0,mlat1 do i=mlon0,mlon1 rhs(i,j) = rhs(i,j)*r0*1.e-2 enddo enddo ! ! rhs here cannot be compared w/ odyn because nhem rhs here is defined ! in nhem indices, whereas in odyn nhem rhs is defined in shem indices. ! However, rhs odyn vs pdyn can be compared after sub gather_pdyn, where ! nhem rhs data is transferred to shem indices of rhs_glb. ! ! call addfld('RHS_N',' ','S',rhs(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) end subroutine rhspde !----------------------------------------------------------------------- subroutine gather_pdyn ! ! Gather needed global arrays to root task, so it can finish non-parallel ! part of dynamo (beginning after sub rhspde) as in original code ! use mpi_module,only: mytid,mp_gather_pdyn ! ! Local: ! 7 fields to gather: zigm11,zigm22,zigmc,zigm2,rim1,rim2,rhs ! integer,parameter :: nf = 7 real :: fmsub(mlon0:mlon1,mlat0:mlat1,nf) real :: fmglb(nmlonp1,nmlat,nf) real :: rhs_nhem(nmlonp1,nmlat) integer :: i,j,jj ! ! These calls exclude halo points in zigm11, etc. ! fmsub(:,:,1) = zigm11 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,2) = zigm22 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,3) = zigmc (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,4) = zigm2 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,5) = rim1 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1) fmsub(:,:,7) = rhs (mlon0:mlon1,mlat0:mlat1) call mp_gather_pdyn(fmsub,mlon0,mlon1,mlat0,mlat1, | fmglb,nmlonp1,nmlat,nf) ! ! Gather 3d zpot in separate call (tell mp_gather_pdyn there ! are nmlevp1) fields. When building under intel/ifort, must have ! -heap-arrays flag set for this call to work: ! if (mytid==0) then zigm11_glb(:,:) = fmglb(:,:,1) zigm22_glb(:,:) = fmglb(:,:,2) zigmc_glb(:,:) = fmglb(:,:,3) zigm2_glb(:,:) = fmglb(:,:,4) rim_glb(:,:,1) = fmglb(:,:,5) rim_glb(:,:,2) = fmglb(:,:,6) rhs_nhem(:,:) = fmglb(:,:,7) ! ! Transfer local global rhs_nhem (from sub rhspde) to rhs_glb, ! so the latter has data in south hemisphere. rhs_glb= 0. ! init do j=1,nmlat ! ! Transfer north pole to equator: if (j == nmlat) then do i=1,nmlonp1 rhs_glb(i,jlateq) = rhs_nhem(i,j) enddo ! Transfer equator to south pole: elseif (j == jlateq) then do i=1,nmlonp1 rhs_glb(i,1) = rhs_nhem(i,j) enddo ! Transfer north hem to south hem: elseif (j > jlateq) then ! 50 -> 96 jj = j-jlateq+1 ! 2 -> 48 do i=1,nmlonp1 rhs_glb(i,jj) = rhs_nhem(i,j) enddo endif enddo ! j=mlat0,mlat1 ! ! This rhs can be compared w/ odyn rhs after rhspde: ! call addfld('RHS',' ','S',rhs_glb,'mlon',1,nmlonp1, ! | 'mlat',1,nmlat,0) endif ! mytid==0 end subroutine gather_pdyn !----------------------------------------------------------------------- subroutine productAX(cofum,phi,mlon0,mlon1,mlat0,mlat1,ax) ! ! am 9/14/2010 ! calculates Ax and put it into ax ! implicit none ! ! Args: integer, intent(in):: mlon0,mlon1,mlat0,mlat1 real, intent(in):: cofum(mlon0:mlon1,mlat0:mlat1,9), | phi(mlon0-1:mlon1+1,mlat0-1:mlat1+1) ! when called with phimsolv ! ! Local: real, intent(out)::ax(mlon0:mlon1,mlat0:mlat1) integer :: i,j,lonend ! ! 9/21/11 btf: Set lonend to avoid exceeding upper bound of x(i+1,..) ! (note eastern most tasks have mlon1==nmlonp1) ! lonend = mlon1 if (mlon1 == nmlon+1) lonend = mlon1-1 do j = mlat0,mlat1 do i = mlon0,lonend ax(i,j) = ( + cofum(i,j,1)*phi(i+1,j)+ + cofum(i,j,2)*phi(i+1,j+1)+ + cofum(i,j,3)*phi(i,j+1)+ + cofum(i,j,4)*phi(i-1,j+1)+ + cofum(i,j,5)*phi(i-1,j)+ + cofum(i,j,6)*phi(i-1,j-1)+ + cofum(i,j,7)*phi(i,j-1)+ + cofum(i,j,8)*phi(i+1,j-1)+ + cofum(i,j,9)*phi(i,j)) enddo enddo ! ! Save product to secondary histories: ! call addfld('AX','AX product',' ',ax,'mlon',mlon0,mlon1, ! | 'mlat',mlat0,mlat1,0) end subroutine productAX !----------------------------------------------------------------------- subroutine highlat_poten ! ! Global PDE solution rim_glb(:,:,1) has been scattered to mag subdomains ! in rim1, and halos set (this overwrites previous rim1 from fieldline ! integrations). Now add high latitude potential from empirical model ! (heelis or weimer), defining phim2d on mag subdomains. After this pthreed ! will expand phim2d to phim3d. ! ! Input: rim1(mag subdomains) ! Solution from mudpack solver (nhem only) ! pfrac(nmlonp1,nmlat0) ! NH fraction of potential ! phihm(nmlonp1,nmlat) ! potential in magnetic ! Output: phim2d(mag subdomains) ! solution with phihm in both nhem and shem ! Both rim1 and phim2d are dimensioned (mlon00:mlon11,mlat00:mlat11) ! ! Both phihm and pfrac have been set by either heelis or weimer. ! phihm is on 2d global mag grid, pfrac is in north hemisphere only ! use input_module,only: iamie use amie_module,only: pcp_nh_amie,pcp_sh_amie use hist_module,only: modeltime use init_module,only: iyear use cons_module,only: dlatm,dlonm,pi_dyn,ylatm,rtd,crit ! Local: logical,parameter :: mod_heelis = .false. ! true == modified integer :: i,j,jn,js real :: fac ! ! Add empirical model potential at high latitude: fac = 1.0 if (mod_heelis) fac = 0. ! modified heelis do j=mlat0,mlat1 if (j > nmlath) cycle ! south only (including equator) jn = nmlat-j+1 js = nmlath-j+1 do i=mlon0,mlon1 phim2d(i,j) = rim1(i,j)+fac*(1.-pfrac(i,js))*(phihm(i,j)- | phihm(i,jn)) enddo enddo do j=mlat0,mlat1 if (j <= nmlath) cycle ! north only (excluding equator) do i=mlon0,mlon1 phim2d(i,j) = rim1(i,j) enddo enddo ! ! There is also a PHIM2D saved at beginning of pthreed below. ! call addfld('PHIM2D',' ','VOLTS', ! | phim2d(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('RIM1',' ','VOLTS', ! | rim1(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) end subroutine highlat_poten !----------------------------------------------------------------------- subroutine pthreed ! ! phim2d is now 2d electric potential on mag subdomains, ! with high-latitude potential added from empirical model (see sub ! highlat_poten), and mag halos set. Now expand phim2d in vertical, ! defining phim3d. Also calculate electric field ed13d, ed23d for ! later current calculations, and ephi3d, elam3d and emz3d for conversion ! to geographic grid (sub pefield), and calculation of ion drifts ! by sub ionvel. ! use cons_module,only: ylatm,h0,re,pi_dyn,table,dlatm,dlonm,rcos0s, | r0,dt1dts,dt0dts use mpi_module,only: mp_mageq_jpm3,mp_magpole_2d,mp_mag_jslot, | mp_mag_halos,mp_magpoles use diags_module,only: mkdiag_ED12 ! ! Local: real,parameter :: eps = 1.e-10, unitvm(nmlon)=1. integer,parameter :: mxneed=nmlat+2 integer :: i,j,k,n,mlon00,mlon11,mlat00,mlat11 real :: csth0,cosltm,sinltm,sym,pi,phims,phimn,siniq real,dimension(nmlonp1) :: thetam,pslot,qslot integer,dimension(nmlonp1) :: islot,jslot,ip1f,ip2f,ip3f real,dimension(mlon0-1:mlon1+1,mlat0-1:mlat1+1) :: ed1,ed2 real,dimension(mlon0-1:mlon1+1,mlat0-1:mlat1+1) :: ephi,elam real :: fpole2d_jpm2(nmlonp1,4,4) ! global lons at S pole+1,2 and N pole-1,2 real :: fpoles(nmlonp1,2,1) ! global lons at poles (1 field only) real :: fmsub(mlon0:mlon1,mlat0:mlat1,4) real :: fmsub1(mlon0-1:mlon1+1,mlat0-1:mlat1+1,5) real :: feq_jpm3(nmlonp1,-3:3,1) ! global lons at equator +/- 3 integer :: jneed(mxneed) ! lats needed from other tasks for interp integer :: njneed,icount real,dimension(mlon0-1:mlon1+1,mxneed) :: | phineed, ! phim2d at needed latitudes | ed1need, ! ed1 at needed latitudes | ed2need, ! ed2 at needed latitudes | ephineed, ! ephi at needed latitudes | elamneed ! elam at needed latitudes real,dimension(mlon0-1:mlon1+1,mxneed,5) :: fmneed real :: phi0j0,phi1j0,phi0j1,phi1j1 real :: ed1i0j0,ed1i1j0,ed1i0j1,ed1i1j1 real :: ed2i0j0,ed2i1j0,ed2i0j1,ed2i1j1 real :: ephi0j0,ephi1j0,ephi0j1,ephi1j1 real :: elam0j0,elam1j0,elam0j1,elam1j1 ! ! External: integer,external :: ixfind real,external :: sddot if (debug) | write(6,"('Enter pthreed: phim2d=',2e12.4)") | minval(phim2d),maxval(phim2d) pi = pi_dyn mlon00=mlon0-1 ; mlon11=mlon1+1 mlat00=mlat0-1 ; mlat11=mlat1+1 ! call addfld('PHIM2D',' ','VOLTS', ! | phim2d(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! ! Calculate ed1,ed2 components of electric field: ! phim2d has halos set, so when mlon0==1, i-1 should wrap ! to value at i==nmlon, and when mlon1==nmlonp1, i+1 should ! wrap to value at i==2. ! ed1 = 0. ephi = 0. do j=mlat0,mlat1 if (j==1.or.j==nmlat) cycle csth0 = cos(-pi/2.+(j-1)*dlatm) do i=mlon0,mlon1 ed1(i,j) = -(phim2d(i+1,j)-phim2d(i-1,j))/(2.*dlonm*csth0)* | rcos0s(j)/(r0*1.e-2) ephi(i,j) = ed1(i,j)*(r0*1.e-2) enddo ! i=mlon0,mlon1 enddo ! j=mlat0,mlat1 ! ! Southern hemisphere (excluding equator): ed2 = 0. elam = 0. do j=mlat0,mlat1 if (j >= nmlath) cycle do i=mlon0,mlon1 ed2(i,j)= -(phim2d(i,j+1)-phim2d(i,j-1))/(2.*dlatm)*dt1dts(j)/ | (r0*1.e-2) elam(i,j)= -(phim2d(i,j+1)-phim2d(i,j-1))/(2.*dlatm)*dt0dts(j) enddo ! i=mlon0,mlon1 enddo ! j=mlat0,mlat1 if (debug) | write(6,"('pthreed after shem: ed1=',2e12.4,' ed2=',2e12.4, | ' ephi=',2e12.4,' elam=',2e12.4)") minval(ed1),maxval(ed1), | minval(ed2),maxval(ed2),minval(ephi),maxval(ephi), | minval(elam),maxval(elam) ! ! Northern hemisphere (excluding equator): do j=mlat0,mlat1 if (j <= nmlath) cycle do i=mlon0,mlon1 ed2(i,j) = (phim2d(i,j+1)-phim2d(i,j-1))/(2.*dlatm)*dt1dts(j)/ | (r0*1.e-2) elam(i,j)= -(phim2d(i,j+1)-phim2d(i,j-1))/(2.*dlatm)*dt0dts(j) enddo ! i=mlon0,mlon1 enddo ! j=mlat0,mlat1 ! call addfld('ED1_a',' ',' ',ed1(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ED2_a',' ',' ',ed2(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('EPHI_a',' ',' ',ephi(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ELAM_a',' ',' ',elam(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! ! Need ed1,2 at global longitudes at j==2 and j==nmlat-1: ! mp_magpole_2d: Return fpole_jpm2(nglblon,1->4,nf) as: ! 1: j = jspole+1 (spole+1) ! 2: j = jspole+2 (spole+2) (unused here) ! 3: j = jnpole-1 (npole-1) ! 4: j = jnpole-2 (npole-2) (unused here) ! fmsub(:,:,1) = ed1(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,2) = ed2(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,3) = ephi(mlon0:mlon1,mlat0:mlat1) fmsub(:,:,4) = elam(mlon0:mlon1,mlat0:mlat1) ! Returns fpole2d_jpm2(nmlonp1,4,4) ! global lons at S pole+1,2 and N pole-1,2 call mp_magpole_2d(fmsub,mlon0,mlon1,mlat0,mlat1,nmlonp1, | 1,nmlat,fpole2d_jpm2,4) ! ! Poles: average over four surrounding points do i = 1,nmlonp1 ip1f(i) = i + nmlon/4 if (ip1f(i) > nmlonp1) ip1f(i) = ip1f(i) - nmlon ip2f(i) = i + nmlon/2 if (ip2f(i) > nmlonp1) ip2f(i) = ip2f(i) - nmlon ip3f(i) = i + 3*nmlon/4 if (ip3f(i) > nmlonp1) ip3f(i) = ip3f(i) - nmlon enddo ! i=1,nmlonp1 ! ! S pole: do j=mlat0,mlat1 if (j==1) then do i=mlon0,mlon1 ed1(i,j)=.25*(fpole2d_jpm2(i,1,1)-fpole2d_jpm2(ip2f(i),1,1)+ | fpole2d_jpm2(ip1f(i),1,2)-fpole2d_jpm2(ip3f(i),1,2)) ed2(i,j)=.25*(fpole2d_jpm2(i,1,2)-fpole2d_jpm2(ip2f(i),1,2)- | fpole2d_jpm2(ip1f(i),1,1)+fpole2d_jpm2(ip3f(i),1,1)) ephi(i,j)=.25*(fpole2d_jpm2(i,1,3)-fpole2d_jpm2(ip2f(i),1,3)+ | fpole2d_jpm2(ip1f(i),1,4)-fpole2d_jpm2(ip3f(i),1,4)) elam(i,j)=.25*(fpole2d_jpm2(i,1,4)-fpole2d_jpm2(ip2f(i),1,4)- | fpole2d_jpm2(ip1f(i),1,3)+fpole2d_jpm2(ip3f(i),1,3)) enddo ! i=mlon0,mlon1 ! ! N pole: elseif (j==nmlat) then do i=mlon0,mlon1 ed1(i,j)=.25*(fpole2d_jpm2(i,3,1)-fpole2d_jpm2(ip2f(i),3,1)+ | fpole2d_jpm2(ip1f(i),3,2)-fpole2d_jpm2(ip3f(i),3,2)) ed2(i,j)=.25*(fpole2d_jpm2(i,3,2)-fpole2d_jpm2(ip2f(i),3,2)- | fpole2d_jpm2(ip1f(i),3,1)+fpole2d_jpm2(ip3f(i),3,1)) ephi(i,j)=.25*(fpole2d_jpm2(i,3,3)-fpole2d_jpm2(ip2f(i),3,3)+ | fpole2d_jpm2(ip1f(i),3,4)-fpole2d_jpm2(ip3f(i),3,4)) elam(i,j)=.25*(fpole2d_jpm2(i,3,4)-fpole2d_jpm2(ip2f(i),3,4)- | fpole2d_jpm2(ip1f(i),3,3)+fpole2d_jpm2(ip3f(i),3,3)) enddo ! i=mlon0,mlon1 endif ! S or N pole enddo ! j=mlat0,mlat1 ! call addfld('ED1_b',' ',' ',ed1(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ED2_b',' ',' ',ed2(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('EPHI_b',' ',' ',ephi(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ELAM_b',' ',' ',elam(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! ! Equator: derivative of quadratic polynomial (3 given points) ! For equator and equator +/- 1 of ed2, we need equator and ! equator +/- 3 of phim2d (note feq_jpm3(nmlonp1,-3:3,1)): ! fmsub(:,:,1) = phim2d(mlon0:mlon1,mlat0:mlat1) call mp_mageq_jpm3(fmsub(:,:,1),mlon0,mlon1,mlat0,mlat1,nmlonp1, | feq_jpm3,1) ! ! Note timegcm does not use dt1dts here, but tiegcm does. do j=mlat0,mlat1 if (j==nmlath-1) then ! equator-1 do i=mlon0,mlon1 ed2(i,j) = (4.*feq_jpm3(i,-2,1)-feq_jpm3(i,-3,1)- | 3.*feq_jpm3(i,-1,1))/(2.*dlatm)*dt1dts(nmlath-1)/(r0*1.e-2) enddo elseif (j==nmlath) then ! equator do i=mlon0,mlon1 ed2(i,j) = (4.*feq_jpm3(i,1,1)-feq_jpm3(i,2,1)- | 3.*feq_jpm3(i,0,1))/(2.*dlatm)*dt1dts(nmlath)/(r0*1.e-2) elam(i,j) = (4.*feq_jpm3(i,1,1)-feq_jpm3(i,2,1)- | 3.*feq_jpm3(i,0,1))/(2.*dlatm) enddo elseif (j==nmlath+1) then ! equator+1 do i=mlon0,mlon1 ed2(i,j) = (4.*feq_jpm3(i,2,1)-feq_jpm3(i,3,1)- | 3.*feq_jpm3(i,1,1))/(2.*dlatm)*dt1dts(nmlath+1)/(r0*1.e-2) enddo endif ! equator +/- 1 enddo ! j=mlat0,mlat1 ! call addfld('ED1_c',' ',' ',ed1(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ED2_c',' ',' ',ed2(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('EPHI_c',' ',' ',ephi(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ELAM_c',' ',' ',elam(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! ! Set halos for 3d calculations: fmsub1(:,:,1) = ed1 fmsub1(:,:,2) = ed2 fmsub1(:,:,3) = ephi fmsub1(:,:,4) = elam if (debug) | write(6,"('pthreed call mp_mag_halos: fmsub1(ed1)=',2e12.4, | ' fmsub1(ed2)=',2e12.4,' fmsub1(ephi)=',2e12.4,' fmsub1(elam)=', | 2e12.4)") minval(fmsub1(:,:,1)),maxval(fmsub1(:,:,1)), | minval(fmsub1(:,:,2)),maxval(fmsub1(:,:,2)), | minval(fmsub1(:,:,3)),maxval(fmsub1(:,:,3)), | minval(fmsub1(:,:,4)),maxval(fmsub1(:,:,4)) call mp_mag_halos(fmsub1,mlon0,mlon1,mlat0,mlat1,1,4) ed1 = fmsub1(:,:,1) ed2 = fmsub1(:,:,2) ephi = fmsub1(:,:,3) elam = fmsub1(:,:,4) ! call addfld('ED1',' ',' ',ed1(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ED2',' ',' ',ed2(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('EPHI',' ',' ',ephi(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call addfld('ELAM',' ',' ',elam(mlon0:mlon1,mlat0:mlat1), ! | 'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! ! Determine latitudes needed for interpolation that fall ! outside a task's latitudinal subdomain: njneed = 0 ! number of unique latitudes needed jneed(:) = -1 ! j-indices of needed latitudes do k=mlev0,mlev1 do j=mlat0,mlat1 if (j==1.or.j==nmlat) cycle ! exclude poles sym = 1. if (j < nmlath) sym = -1. cosltm = cos(ylatm(j)) do i=mlon0,mlon1 if (i==nmlonp1) cycle thetam(i)=(re+zpotm3d(i,j,mlev0))/(re+zpotm3d(i,j,k)) ! write(6,"('pthreed: i=',i4,' j=',i4,' k=',i4, ! | ' zpotm3d(i,j,mlev0)=',e12.4,' zpotm3d(i,j,k)=',e12.4, ! | ' thetam(i)=',e12.4)") i,j,k,zpotm3d(i,j,mlev0), ! | zpotm3d(i,j,k),thetam(i) thetam(i) = acos(sqrt(thetam(i))*cosltm*(1.-eps)) ! write(6,"('pthreed: i=',i4,' j=',i4,' k=',i4,' cosltm=', ! | e12.4,' thetam(i)=',e12.4)") i,j,k,cosltm,thetam(i) pslot(i) = thetam(i)*180./pi+1. islot(i) = pslot(i) pslot(i) = pslot(i)-float(islot(i)) thetam(i) = ((1.-pslot(i))*table(islot(i),2)+pslot(i)* | table(islot(i)+1,2))*sym ! thetam negative for south hem islot(i) = i pslot(i) = 0. qslot(i) = (thetam(i)+pi/2.)/dlatm+1. jslot(i) = qslot(i) qslot(i) = qslot(i)-float(jslot(i)) ! ! Save j index if outside subdomain w/ halos: if ((jslot(i) < mlat00 .or. jslot(i) > mlat11).and. | .not.any(jslot(i)==jneed)) then njneed = njneed+1 if (njneed > mxneed) call shutdown('njneed') jneed(njneed) = jslot(i) endif ! jslot is outside subdomain ! ! Save j+1 index if outside subdomain: if ((jslot(i)+1 < mlat00 .or. jslot(i)+1 > mlat11).and. | .not.any(jslot(i)+1==jneed)) then njneed = njneed+1 if (njneed > mxneed) call shutdown('njneed') jneed(njneed) = jslot(i)+1 endif ! jslot(i)+1 is outside subdomain enddo ! i=mlon0,mlon1 enddo ! j=mlat0,mlat1 enddo ! k=mlev0,mlev1 ! ! Get phim2 at needed latitudes (note inclusion of phim2d halos). ! real,intent(in) :: fin(mlon00:mlon11,mlat00:mlat11,nf) ! data at current subdomain ! real,intent(out) :: fout(mlon00:mlon11,mxneed,nf) ! returned data at needed lats ! fmsub1(:,:,1) = phim2d fmsub1(:,:,2) = ed1 fmsub1(:,:,3) = ed2 fmsub1(:,:,4) = ephi fmsub1(:,:,5) = elam if (debug) write(6,"('pthreed call mp_mag_jslot')") call mp_mag_jslot(fmsub1,mlon00,mlon11,mlat00,mlat11, | fmneed,jneed,mxneed,5) if (debug) write(6,"('pthreed after mp_mag_jslot')") phineed = fmneed(:,:,1) ed1need = fmneed(:,:,2) ed2need = fmneed(:,:,3) ephineed= fmneed(:,:,4) elamneed= fmneed(:,:,5) ! ephi3d = 0. elam3d = 0. emz3d = 0. do k=mlev0,mlev1 do j=mlat0,mlat1 if (j==1.or.j==nmlat) cycle ! exclude poles sym = 1. if (j < nmlath) sym = -1. cosltm = cos(ylatm(j)) sinltm = sin(ylatm(j)) do i=mlon0,mlon1 if (i==nmlonp1) cycle siniq = 2.*sinltm/sqrt(4.-3.*cosltm**2) thetam(i)=(re+zpotm3d(i,j,mlev0))/(re+zpotm3d(i,j,k)) thetam(i) = acos(sqrt(thetam(i))*cosltm*(1.-eps)) pslot(i) = thetam(i)*180./pi+1. islot(i) = pslot(i) pslot(i) = pslot(i)-float(islot(i)) thetam(i) = ((1.-pslot(i))*table(islot(i),2)+pslot(i)* | table(islot(i)+1,2))*sym ! thetam negative for south hem islot(i) = i pslot(i) = 0. qslot(i) = (thetam(i)+pi/2.)/dlatm+1. jslot(i) = qslot(i) qslot(i) = qslot(i)-float(jslot(i)) ! ! Check for jslot in subdomain: if (jslot(i) >= mlat00.and.jslot(i) <= mlat11) then ! within subdomain phi0j0 = phim2d(islot(i) ,jslot(i)) phi1j0 = phim2d(islot(i)+1,jslot(i)) ed1i0j0 = ed1(islot(i) ,jslot(i)) ed1i1j0 = ed1(islot(i)+1,jslot(i)) ed2i0j0 = ed2(islot(i) ,jslot(i)) ed2i1j0 = ed2(islot(i)+1,jslot(i)) ephi0j0 = ephi(islot(i) ,jslot(i)) ephi1j0 = ephi(islot(i)+1,jslot(i)) elam0j0 = elam(islot(i) ,jslot(i)) elam1j0 = elam(islot(i)+1,jslot(i)) else ! jslot outside subdomain n = ixfind(jneed,mxneed,jslot(i),icount) if (n==0) then write(6,"('>>> pthreed: i=',i4,' j=',i4,' k=',i4, | ' could not find jslot ',i4,' in jneed=',/, | (10i4))") i,j,k,jslot(i),jneed call shutdown('jslot(i) not in jneed') endif phi0j0 = phineed(islot(i) ,n) phi1j0 = phineed(islot(i)+1,n) ed1i0j0 = ed1need(islot(i) ,n) ed1i1j0 = ed1need(islot(i)+1,n) ed2i0j0 = ed2need(islot(i) ,n) ed2i1j0 = ed2need(islot(i)+1,n) ephi0j0 =ephineed(islot(i) ,n) ephi1j0 =ephineed(islot(i)+1,n) elam0j0 =elamneed(islot(i) ,n) elam1j0 =elamneed(islot(i)+1,n) endif ! ! Check for jslot+1 in subdomain: if (jslot(i)+1 >= mlat00.and.jslot(i)+1 <= mlat11) then ! within subdomain phi0j1 = phim2d(islot(i) ,jslot(i)+1) phi1j1 = phim2d(islot(i)+1,jslot(i)+1) ed1i0j1 = ed1(islot(i) ,jslot(i)+1) ed1i1j1 = ed1(islot(i)+1,jslot(i)+1) ed2i0j1 = ed2(islot(i) ,jslot(i)+1) ed2i1j1 = ed2(islot(i)+1,jslot(i)+1) ephi0j1 = ephi(islot(i) ,jslot(i)+1) ephi1j1 = ephi(islot(i)+1,jslot(i)+1) elam0j1 = elam(islot(i) ,jslot(i)+1) elam1j1 = elam(islot(i)+1,jslot(i)+1) else ! jslot+1 outside subdomain n = ixfind(jneed,mxneed,jslot(i)+1,icount) if (n==0) then write(6,"('>>> pthreed: i=',i4,' j=',i4,' k=',i4, | ' could not find jslot+1 ',i4,' in jneed=',/, | (10i4))") i,j,k,jslot(i)+1,jneed call shutdown('jslot(i)+1 not in jneed') endif phi0j1 = phineed(islot(i) ,n) phi1j1 = phineed(islot(i)+1,n) ed1i0j1 = ed1need(islot(i) ,n) ed1i1j1 = ed1need(islot(i)+1,n) ed2i0j1 = ed2need(islot(i) ,n) ed2i1j1 = ed2need(islot(i)+1,n) ephi0j1 =ephineed(islot(i) ,n) ephi1j1 =ephineed(islot(i)+1,n) elam0j1 =elamneed(islot(i) ,n) elam1j1 =elamneed(islot(i)+1,n) endif ! ! Do the interpolation: phim3d(i,j,k) = (1.-qslot(i))*((1.-pslot(i))* | phi0j0 + pslot(i) * phi1j0) + qslot(i)*((1.-pslot(i))* | phi0j1 + pslot(i) * phi1j1) ed13d(i,j,k) = (1.-qslot(i))*((1.-pslot(i))* | ed1i0j0 + pslot(i) * ed1i1j0) + qslot(i)*((1.-pslot(i))* | ed1i0j1 + pslot(i) * ed1i1j1) ed23d(i,j,k) = (1.-qslot(i))*((1.-pslot(i))* | ed2i0j0 + pslot(i) * ed2i1j0) + qslot(i)*((1.-pslot(i))* | ed2i0j1 + pslot(i) * ed2i1j1) ! ephi3d(i,j,k) = (1.-qslot(i))*((1.-pslot(i))* | ephi0j0 + pslot(i) * ephi1j0) + qslot(i)*((1.-pslot(i))* | ephi0j1 + pslot(i) * ephi1j1) elam3d(i,j,k) = (1.-qslot(i))*((1.-pslot(i))* | elam0j0 + pslot(i) * elam1j0) + qslot(i)*((1.-pslot(i))* | elam0j1 + pslot(i) * elam1j1) enddo ! i=mlon0,mlon1 enddo ! j=mlat0,mlat1 enddo ! k=mlev0,mlev1 ! ! Mag poles for phim: ! mp_magpoles returns global longitudes at S,N poles in fpoles(nglblon,2,nf) ! if (debug) write(6,"('pthreed call mp_magpoles')") call mp_magpoles(phim2d(mlon0:mlon1,mlat0:mlat1), | mlon0,mlon1,mlat0,mlat1,nmlonp1,1,nmlat,fpoles,1) if (debug) write(6,"('pthreed after mp_magpoles')") phims=sddot(nmlon,unitvm,fpoles(1:nmlon,1,1))/float(nmlon) phimn=sddot(nmlon,unitvm,fpoles(1:nmlon,2,1))/float(nmlon) do k=mlev0,mlev1 do j=mlat0,mlat1 if (j==1) then do i=mlon0,mlon1 phim3d(i,j,k) = phims ed13d(i,j,k) = ed1(i,j) ed23d(i,j,k) = ed2(i,j) ephi3d(i,j,k) = ephi(i,j) elam3d(i,j,k) = ed2(i,j)*(r0*1.e-2) enddo elseif (j==nmlat) then do i=mlon0,mlon1 phim3d(i,j,k) = phimn ed13d(i,j,k) = ed1(i,j) ed23d(i,j,k) = ed2(i,j) ephi3d(i,j,k) = ephi(i,j) elam3d(i,j,k) = -ed2(i,j)*(r0*1.e-2) enddo endif ! poles enddo ! j=mlat0,mlat1 enddo ! k=mlev0,mlev1 ! ! Upward electric field: do k=mlev0+1,mlev1-1 do j=mlat0,mlat1 do i=mlon0,mlon1 emz3d(i,j,k) = -(phim3d(i,j,k+1)-phim3d(i,j,k-1)) enddo enddo enddo do j=mlat0,mlat1 do i=mlon0,mlon1 emz3d(i,j,mlev1) = emz3d(i,j,mlev1-1) emz3d(i,j,mlev0) = emz3d(i,j,mlev0+1) enddo enddo if (debug) write(6,"('pthreed before mp_mag_periodic_f2d')") do k=mlev0,mlev1 call mp_mag_periodic_f2d(phim3d(:,:,k),mlon0,mlon1, | mlat0,mlat1,1) call mp_mag_periodic_f2d(ed13d(:,:,k),mlon0,mlon1, | mlat0,mlat1,1) call mp_mag_periodic_f2d(ed23d(:,:,k),mlon0,mlon1, | mlat0,mlat1,1) call mp_mag_periodic_f2d(ephi3d(:,:,k),mlon0,mlon1, | mlat0,mlat1,1) call mp_mag_periodic_f2d(ephi3d(:,:,k),mlon0,mlon1, | mlat0,mlat1,1) call mp_mag_periodic_f2d(elam3d(:,:,k),mlon0,mlon1, | mlat0,mlat1,1) call mp_mag_periodic_f2d(emz3d(:,:,k),mlon0,mlon1, | mlat0,mlat1,1) enddo if (debug) write(6,"('pthreed after mp_mag_periodic_f2d')") ! ! Cannot use mlev0 (-2) in addfld call, so plot from level 1: ! do j=mlat0,mlat1 call addfld('PHIM3D',' ',' ',phim3d(:,j,mlev0:mlev1-1), | 'mlon',mlon0,mlon1,'mlev',1,mlevtop,j) call addfld('EMX',' ',' ',ephi3d(:,j,mlev0:mlev1-1), | 'mlon',mlon0,mlon1,'mlev',1,mlevtop,j) call addfld('EMY',' ',' ',elam3d(:,j,mlev0:mlev1-1), | 'mlon',mlon0,mlon1,'mlev',1,mlevtop,j) call addfld('EMZ',' ',' ',emz3d(:,j,mlev0:mlev1-1), | 'mlon',mlon0,mlon1,'mlev',1,mlevtop,j) call mkdiag_ED12('ED1',ed13d(:,j,mlev0:mlev1-1), | mlon0,mlon1,1,mlevtop,j) call mkdiag_ED12('ED2',ed23d(:,j,mlev0:mlev1-1), | mlon0,mlon1,1,mlevtop,j) enddo if (debug) write(6,"('pthreed returning')") end subroutine pthreed !----------------------------------------------------------------------- subroutine pefield ! ! Parallel version of efield. Get electric field on geographic subdomains. ! (ex,ey,ez in this module). This is called from advance, early in the ! timestep. If first timestep, calculate electric field in mag coords from ! phim3d, otherwise use electric field on mag subdomains from pthreed of ! the previous step. ! use cons_module,only: dt0dts,dlatm,dlonm,pi,rcos0s,table use mpi_module,only: mp_magpole_3d,mp_mag_halos,mp_magpoles, | mytid use magfield_module,only: alatm,alonm ! (nlonp1,0:nlatp1) use my_esmf,only: mag_ephi3d,mag_elam3d,mag_emz3d, | geo_ephi3d,geo_elam3d,geo_emz3d use diags_module,only: mkdiag_EXYZ implicit none ! ! Local: integer :: i,ii,j,k,lonbeg,lonend real :: fpole3d_jpm2(nmlonp1,4,mlev0:mlev1,1) ! global lons at S pole+1,2 and N pole-1,2 real :: fpoles(nmlonp1,2,mlev0:mlev1) ! global lons at poles real :: phi3d(mlon0-1:mlon1+1,mlat0-1:mlat1+1,mlev0:mlev1) real :: csth0 integer,dimension(nlonp4) :: islot,jslot real,dimension(nlonp4) :: pslot,qslot,thetam real,dimension(lon0:lon1,lat0:lat1,mlev0:mlev1) :: | exgeo,eygeo,ezgeo if (debug) write(6,"('Enter pefield: istep=',i3)") istep ! do j=mlat0,mlat1 ! call addfld('PHIM3D_EFLD',' ','VOLTS', ! | phim3d(mlon0:mlon1,j,1:mlev1), ! | 'mlon',mlon0,mlon1,'imlev',1,mlev1,j) ! enddo ! ! Note if istep > 1, then 3d electric field ephi3d,elam3d,emz3d was ! calculated by pthreed in previous step. if (istep==1) then if (debug) | write(6,"('pefield (step=1): phim3d(:,:,:) min,max=', | 2e12.4)") minval(phim3d(:,:,:)),maxval(phim3d(:,:,:)) ! ! Copy phim3d to local phi3d, and set halo points: ! (phim3d does not have halo points) ! phi3d = 0. do k=mlev0,mlev1 do j=mlat0,mlat1 do i=mlon0,mlon1 phi3d(i,j,k) = phim3d(i,j,k) enddo enddo enddo do k=mlev0,mlev1 if (debug) write(6,"('pefield (step=1) k=',i4, | ' call mp_mag_halos for phi3d')") k call mp_mag_halos(phi3d(:,:,k),mlon0,mlon1,mlat0,mlat1,1,1) if (debug) write(6,"('pefield (step=1) k=',i4, | ' after mp_mag_halos for phi3d')") k enddo ! ! Return fpole3d_jpm2(nglblon,1->4,mlev0:mlev1,nf) as: ! 1: j = jspole+1 (spole+1) ! 2: j = jspole+2 (spole+2) not used here ! 3: j = jnpole-1 (npole-1) ! 4: j = jnpole-2 (npole-2) not used here ! if (debug) write(6,"('pefield (step=1) call mp_magpole_3d', | ' for phi3d')") call mp_magpole_3d(phim3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1), | mlon0,mlon1,mlat0,mlat1,nmlevp1,nmlonp1,1,nmlat,fpole3d_jpm2,1) if (debug) write(6,"('pefield (step=1) after mp_magpole_3d', | ' for phi3d')") ! ! Set j=0 and j=nmlat+1 of local phi3d. This overwrites the far ! north and south halo points set by mp_mag_halos above. ! do j=mlat0,mlat1 if (j==1) then do i=mlon0,mlon1 ii = 1+mod(i-1+nmlon/2,nmlon) ! over the south pole phi3d(i,j-1,:) = fpole3d_jpm2(ii,1,:,1) enddo elseif (j==nmlat) then do i=mlon0,mlon1 ii = 1+mod(i-1+nmlon/2,nmlon) ! over the north pole phi3d(i,j+1,:) = fpole3d_jpm2(ii,3,:,1) enddo endif ! poles or not enddo ! j=mlat0,mlat1 ! ! Meridional component of electric field: do k=mlev0,mlev1 do j=mlat0,mlat1 do i=mlon0,mlon1 elam3d(i,j,k) = -(phi3d(i,j+1,k)-phi3d(i,j-1,k))/ | (2.*dlatm)*dt0dts(j) enddo enddo enddo ! ! Zonal component of electric field: do j=mlat0,mlat1 if (j==1.or.j==nmlat) cycle csth0 = cos(-pi/2.+float(j-1)*dlatm) do k=mlev0,mlev1 do i=mlon0,mlon1 ephi3d(i,j,k) = -(phi3d(i+1,j,k)-phi3d(i-1,j,k))/ | (2.*dlonm*csth0)*rcos0s(j) enddo enddo enddo ! ! Polar values for ephi3d (need global lons at poles of elam3d): ! if (debug) write(6,"('pefield call mp_magpoles for elam3d')") call mp_magpoles(elam3d(:,:,mlev0:mlev1),mlon0,mlon1, | mlat0,mlat1,nmlonp1,1,nmlat,fpoles,nmlevp1) if (debug) write(6,"('pefield after mp_magpoles for elam3d')") do k=mlev0,mlev1 do j=mlat0,mlat1 if (j==1) then ! south pole do i=mlon0,mlon1 ii = 1+mod(i-1+(nmlon/4),nmlon) ! over the south pole ephi3d(i,j,k) = fpoles(ii,1,k) enddo elseif (j==nmlat) then ! north pole do i=mlon0,mlon1 ii = 1+mod(i-1+((3*nmlon)/4),nmlon) ! over the north pole ephi3d(i,j,k) = fpoles(ii,2,k) enddo endif ! poles or not enddo ! j=mlat0,mlat1 enddo ! k=mlev0,mlev1 ! ! emz = d(phi)/dz do k=mlev0+1,mlev1-1 do j=mlat0,mlat1 do i=mlon0,mlon1 emz3d(i,j,k) = -(phim3d(i,j,k+1)-phi3d(i,j,k-1)) enddo enddo enddo ! k=mlev0+1,mlev1-1 endif ! istep == 1 ! ! If istep > 1, then ephi3d, elam3d and emz3d were set by pthreed ! in the previous step. ! ! do j=mlat0,mlat1 ! call addfld('EMX_EFLD','EMX in pefield',' ',ephi3d(:,j,1:mlev1), ! | 'mlon',mlon0,mlon1,'mlev',1,mlev1,j) ! call addfld('EMY_EFLD','EMY in pefield',' ',elam3d(:,j,1:mlev1), ! | 'mlon',mlon0,mlon1,'mlev',1,mlev1,j) ! call addfld('EMZ_EFLD','EMZ in pefield',' ',emz3d(:,j,1:mlev1), ! | 'mlon',mlon0,mlon1,'mlev',1,mlev1,j) ! enddo ! j=1,nmlat ! ! Use ESMF to regrid the electric field to the geographic grid: ! allocate(ephi3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) ! real,dimension(lon0:lon1,lat0:lat1,mlev0:mlev1) :: exgeo,... if (debug) write(6,"('pefield call mag2geo_3d')") call mag2geo_3d(ephi3d,exgeo,mag_ephi3d,geo_ephi3d,'EPHI3D ') call mag2geo_3d(elam3d,eygeo,mag_elam3d,geo_elam3d,'ELAM3D ') call mag2geo_3d(emz3d, ezgeo,mag_emz3d ,geo_emz3d ,'EMZ3D ') if (debug) write(6,"('pefield after mag2geo_3d')") ! ! Define ex,ey,ez on geographic subdomains for ionvel: do k=1,nlevp1 do j=lat0,lat1 do i=lon0,lon1 ex(k,i,j) = exgeo(i,j,k) ey(k,i,j) = eygeo(i,j,k) ez(k,i,j) = ezgeo(i,j,k) enddo enddo if (debug) | write(6,"('pefield: k=',i4,' ex=',2e12.4,' ey=',2e12.4, | ' ez=',2e12.4)") k, | minval(ex(k,lon0:lon1,lat0:lat1)), | maxval(ex(k,lon0:lon1,lat0:lat1)), | minval(ey(k,lon0:lon1,lat0:lat1)), | maxval(ey(k,lon0:lon1,lat0:lat1)), | minval(ez(k,lon0:lon1,lat0:lat1)), | maxval(ez(k,lon0:lon1,lat0:lat1)) enddo ! k=1,nlevp1 do j=lat0,lat1 call mkdiag_EXYZ('EX',ex(1:nlevp1,lon0:lon1,j),1,nlevp1, | lon0,lon1,j) call mkdiag_EXYZ('EY',ey(1:nlevp1,lon0:lon1,j),1,nlevp1, | lon0,lon1,j) call mkdiag_EXYZ('EZ',ez(1:nlevp1,lon0:lon1,j),1,nlevp1, | lon0,lon1,j) enddo if (debug) write(6,"('pefield returning')") end subroutine pefield !----------------------------------------------------------------------- subroutine mag2geo_3d(fmag,fgeo,ESMF_mag,ESMF_geo,fname) ! From pefield: ! call mag2geo_3d(ephi3d,exgeo,mag_ephi3d,geo_ephi3d,'EPHI3D ') ! allocate(ephi3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat) ! real,dimension(lon0:lon1,lat0:lat1,mlev0:mlev1) :: exgeo,... ! ! Convert field on geomagnetic grid fmag to geographic grid in fgeo. ! use my_esmf,only: esmf_set3d_mag,esmf_regrid,esmf_get_3dfield, | ESMF_Field ! ! Args: character(len=*) :: fname type(ESMF_Field),intent(inout) :: ESMF_mag, ESMF_geo real,intent(in) :: fmag(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1) real,intent(out) :: fgeo(lon0:lon1,lat0:lat1,mlev0:mlev1) ! ! Local: integer :: j character(len=8) :: fnames(1) type(ESMF_Field) :: magfields(1) real,pointer,dimension(:,:,:) :: fptr fgeo = 0. fnames(1) = fname magfields(1) = ESMF_mag ! ! Put fmag into ESMF mag field on mag source grid: call esmf_set3d_mag(magfields,fnames,fmag,1,mlev0,mlev1, | mlon0,mlon1,mlat0,mlat1) ! ! Regrid to geographic destination grid, defining ESMF_geo: call esmf_regrid(ESMF_mag,ESMF_geo,'mag2geo',3) ! ! Put regridded geo field into pointer: call esmf_get_3dfield(ESMF_geo,fptr,fname) ! ! Transfer from pointer to output arg: do j=lat0,lat1 fgeo(:,j,:) = fptr(:,j,:) enddo end subroutine mag2geo_3d !----------------------------------------------------------------------- ! #if defined(INTERCOMM) || defined(CISMAH) subroutine mag2geo_2d(fmag,fgeo,ESMF_mag,ESMF_geo,fname) ! ! Convert field on geomagnetic grid fmag to geographic grid in fgeo. ! use my_esmf,only: esmf_set2d_mag,esmf_regrid,esmf_get_2dfield, | ESMF_Field ! ! Args: character(len=*) :: fname type(ESMF_Field),intent(inout) :: ESMF_mag, ESMF_geo real,intent(in) :: fmag(mlon0:mlon1,mlat0:mlat1) real,intent(out) :: fgeo(lon0:lon1,lat0:lat1) ! ! Local: integer :: j character(len=8) :: fnames(1) type(ESMF_Field) :: magfields(1) real,pointer,dimension(:,:) :: fptr fgeo = 0. fnames(1) = fname magfields(1) = ESMF_mag ! ! Put fmag into ESMF mag field on mag source grid: call esmf_set2d_mag(magfields,fnames,fmag,1, | mlon0,mlon1,mlat0,mlat1) ! ! Regrid to geographic destination grid, defining ESMF_geo: call esmf_regrid(ESMF_mag,ESMF_geo,'mag2geo',2) ! ! Put regridded geo field into pointer: call esmf_get_2dfield(ESMF_geo,fptr,fname) ! ! Transfer from pointer to output arg: do j=lat0,lat1 fgeo(:,j) = fptr(:,j) enddo end subroutine mag2geo_2d #endif !----------------------------------------------------------------------- subroutine geo2mag_3d(fgeo,fmag,ESMF_geo,ESMF_mag,fname) ! ! Convert field on geographic grid fgeo to geomagnetic grid in fmag. ! use my_esmf,only: esmf_set3d_geo,esmf_regrid,esmf_get_3dfield, | ESMF_Field ! ! Args: character(len=*) :: fname type(ESMF_Field),intent(inout) :: ESMF_geo, ESMF_mag real,intent(in) :: fgeo(nlevp1,lon0:lon1,lat0:lat1) real,intent(out) :: fmag(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1) ! ! Local: integer :: j,k character(len=8) :: fnames(1) type(ESMF_Field) :: geofields(1) real,pointer,dimension(:,:,:) :: fptr real :: fgeotmp(mlev0:mlev1,lon0:lon1,lat0:lat1) integer :: lbnd(3),ubnd(3) ! do j=lat0,lat1 ! write(6,"('geo2mag_3d: j=',i4,' fgeo(:,lon0:lon1,j)=', ! | 2e12.4)") j,minval(fgeo(:,lon0:lon1,j)), ! | maxval(fgeo(:,lon0:lon1,j)) ! enddo fmag = 0. fnames(1) = fname geofields(1) = ESMF_geo ! ! Put fgeo into ESMF geo field on geo source grid: ! subroutine esmf_set3d_geo(fields,fnames,f,nf,ilev0,ilev1, ! ilon0,ilon1,ilat0,ilat1) ! ! This leaves bottom k-levels (-2 to 0 or -6 to 0 depending on resolution) ! of fgeotmp zeroed out: ! fgeotmp = 0. fgeotmp(1:nlevp1,:,:) = fgeo(:,:,:) ! ! Note esmf_set3d_geo will change dimension order from input ! fgeotmp(k,i,j) to output ESMF_geo(i,j,k). ! call esmf_set3d_geo(geofields,fnames,fgeotmp,1,mlev0,mlev1, | lon0,lon1,lat0,lat1) ! ! Get ESMF_geo and print bounds and min,max of levels 1,nlevp1: call esmf_get_3dfield(ESMF_geo,fptr,fname) ! write(6,"('geo2mag_3d after esmf_get_3dfield: lbnd(fptr)=',3i4, ! | ' ubnd(fptr)=',3i4)") lbound(fptr),ubound(fptr) ! do j=lat0,lat1 ! write(6,"('geo2mag_3d: j=',i4,' fptr to ', ! | 'ESMF_geo(:,j,1:nlevp1)=',2e12.4)") ! | j,minval(fptr(:,j,1:nlevp1)),maxval(fptr(:,j,1:nlevp1)) ! enddo ! ! Regrid to magnetic destination grid, defining ESMF_mag: call esmf_regrid(ESMF_geo,ESMF_mag,'geo2mag',3) ! ! Put regridded mag field into pointer: call esmf_get_3dfield(ESMF_mag,fptr,fname) ! ! Transfer from pointer to output arg: do j=mlat0,mlat1 fmag(:,j,:) = fptr(:,j,:) ! write(6,"('geo2mag_3d: j=',i4,' fmag(:,j,1:mlev1)=',2e12.4)") ! | j,minval(fmag(:,j,1:mlev1)),maxval(fmag(:,j,1:mlev1)) enddo end subroutine geo2mag_3d !----------------------------------------------------------------------- subroutine prepare_phig3d(wrprim,wrsech) ! ! If its time to write a history, convert phim3d to geographic grid in ! phig3d, and transfer to appropriate f4d array, where it will be ! gathered to the root task before writing to the history. ! use fields_module,only: itc,poten ! pointer to f4d(i_poten) use my_esmf,only: mag_phi3d,geo_phi3d ! ! Args: logical,intent(in) :: wrprim,wrsech ! ! Local: integer :: i,j ! if (.not.wrprim.and..not.wrsech) return ! ! Transform 3d mag phim3d to geo phig3d, for write to histories: call mag2geo_3d(phim3d,phig3d,mag_phi3d,geo_phi3d,'PHIM3D ') ! ! do j=lat0,lat1 ! call addfld('PHIG3D','3D EPOT on Geographic grid','Volts', ! | phig3d(lon0:lon1,j,:),'lon',lon0,lon1,'ilev',1,nlevp1,j) ! enddo ! j=lat0,lat1 ! ! Update poten (pointer to f4d(i_poten)) for gather2root and write to ! history. In fields.F: poten(levd0:levd1,lond0:lond1,latd0:latd1,2) ! do j=lat0,lat1 do i=lon0,lon1 poten(:,i,j,itc) = phig3d(i,j,:) enddo enddo ! do j=lat0,lat1 ! call addfld('POTEN_PHIG3D','3D EPOT on Geographic grid','Volts', ! | poten(:,lon0:lon1,j,itc),'ilev',1,nlevp1,'lon',lon0,lon1,j) ! enddo ! j=lat0,lat1 end subroutine prepare_phig3d !----------------------------------------------------------------------- subroutine define_phim3d(ixt) ! ! Convert electric potential read from source history on geographic ! grid to magnetic grid in phim3d. This is called once per run from ! sub readsource (rdsource.F). ! use my_esmf,only: mag_ped,geo_ped use fields_module,only: poten ! ! Args: integer,intent(in) :: ixt ! ! Local: integer :: j ! ! Define phim3d from poten if doing dynamo: ! (poten was read from source history) ! do j=lat0,lat1 ! write(6,"('define_phim3d: j=',i4,' poten(:,lon0:lon1,j,ixt)=', ! | 2e12.4)") j,minval(poten(:,lon0:lon1,j,ixt)), ! | maxval(poten(:,lon0:lon1,j,ixt)) ! enddo call geo2mag_3d(poten(:,lon0:lon1,lat0:lat1,ixt), | phim3d,geo_ped,mag_ped,'PHIG3D ') ! do j=mlat0,mlat1 ! write(6,"('define_phim3d: j=',i3,' phim3d(:,j,:) min,max=', ! | 2e12.4)") j,minval(phim3d(:,j,:)),maxval(phim3d(:,j,:)) ! enddo end subroutine define_phim3d !----------------------------------------------------------------------- ! ! May, 2012 btf: ! From here to end of module are routines from original serial code, ! which operate in global magnetic space, i.e., the remaining serial ! portion of the dynamo code. These routine will either be replaced by ! a new parallel solver, or they will be parallelized in shared memory ! with OpenMP directives. ! !----------------------------------------------------------------------- subroutine stencils use cons_module,only: pi_dyn,dlatm,dlonm ! ! Locals: integer :: i,j,jj,jjj,j0,n,ncc,nmaglon,nmaglat real :: sym real :: cs(nmlat0) real :: c0_plt(nmlonp1,nmlat,10) real :: cofum_plt(nmlonp1,nmlat,9) ! ! Set index array nc and magnetic latitude cosine array: ! nc pointes to the start of the coefficient array for each level nc(1) = nc0 nc(2) = nc1 nc(3) = nc2 nc(4) = nc3 nc(5) = nc4 nc(6) = ncee do j=1,nmlat0 cs(j) = cos(pi_dyn/2.-(nmlat0-j)*dlatm) enddo ! j=1,nmlat0 ! ! Set up difference coefficients. Replace zigm11 by A, zigm22 by B, ! zigmc by C, and zigm2 by D. ! j0 = nmlat0-nmlath do j=1,nmlath ! 1,49 (assuming nmlat=97) jj = nmlath+j-1 ! 49,97 jjj = nmlath-j+1 ! 49,1 ! ! factor 4 from 5-point diff. stencil ! Sigma_(phi lam)/( 4*Delta lam* Delta lon ) ! Sigma_(phi lam)/( 4*Delta lam* Delta lon ) ! Sigma_(lam lam)*cos(lam_m)*DT0DTS/(Delta lam)^2 ! -zigmc_north = southern hemis. 49,1 equator-pole ! -zigm2_north = southern hemis. 49,1 equator-pole ! zigm22 = southern hemis. 49,1 equator-pole ! do i=1,nmlonp1 zigmc_glb(i,jj) = (zigmc_glb(i,jj) +zigm2_glb(i,jj))/ | (4.*dlatm*dlonm) zigm2_glb(i,jj) = zigmc_glb(i,jj)-2.*zigm2_glb(i,jj)/ | (4.*dlatm*dlonm) zigm22_glb(i,jj) = zigm22_glb(i,jj)*cs(j0+j)/dlatm**2 zigmc_glb(i,jjj) = -zigmc_glb(i,jj) zigm2_glb(i,jjj) = -zigm2_glb(i,jj) zigm22_glb(i,jjj) = zigm22_glb(i,jj) enddo ! i=1,nmlonp1 if (j /= nmlath) then ! ! Sigma_(phi phi)/( cos(lam_m)*DT0DTS*(Delta lon)^2 ) ! zigm11 = southern hemis. 49,1 equator-pole ! do i = 1,nmlonp1 zigm11_glb(i,jj) = zigm11_glb(i,jj)/(cs(j0+j)*dlonm**2) zigm11_glb(i,jjj) = zigm11_glb(i,jj) enddo endif enddo ! j=1,nmlath ! ! Set zigm11 to zero at megnetic poles to avoid floating exception ! (values at poles are not used): ! do i = 1,nmlonp1 zigm11_glb(i,1) = 0.0 zigm11_glb(i,nmlat) = 0.0 enddo ! ! Save global zigm for comparison w/ odyn: ! real,dimension(nmlonp1,nmlat) :: zigm11_glb,zigm22_glb,zigmc_glb,zigm2_glb,rhs_glb ! ! call addfld('ZIGM11_GLB',' ',' ',zigm11_glb,'mlon',1,nmlonp1, ! | 'mlat',1,nmlat,0) ! call addfld('ZIGM22_GLB',' ',' ',zigm22_glb,'mlon',1,nmlonp1, ! | 'mlat',1,nmlat,0) ! call addfld('ZIGMC_GLB',' ',' ',zigmc_glb,'mlon',1,nmlonp1, ! | 'mlat',1,nmlat,0) ! call addfld('ZIGM2_GLB',' ',' ',zigm2_glb,'mlon',1,nmlonp1, ! | 'mlat',1,nmlat,0) ! call addfld('RHS_GLB',' ',' ',rhs_glb,'mlon',1,nmlonp1, ! | 'mlat',1,nmlat,0) ! ! Clear array for difference stencils at all levels: call clearcee(cee,nmlon0,nmlat0) ! ! isolve = 2 -> modified mudpack solver (modified and unmodified coefficients) ! (currently the only option available) ! ! real,dimension(nmlon0,nmlat0,9) :: cofum cofum(:,:,:) = 0. ! init ! ! Calculate contribution to stencils from each PDE coefficient ! ! Sigma_(phi phi)/( cos(lam_m)*dt0dts*(Delta lon)^2 ) sym = 1. call stencmd(zigm11_glb,cs,nmlon0,nmlat0,sym,cee,1) ! ! Sigma_(lam lam)*cos(lam_m)*dt0dts/(Delta lam)^2 sym = 1. call stencmd(zigm22_glb,cs,nmlon0,nmlat0,sym,cee,4) ! ! Sigma_(phi lam)/( 4*Delta lam* Delta lon ) sym = -1. call stencmd(zigmc_glb,cs,nmlon0,nmlat0,sym,cee,2) ! ! Sigma_(lam phi)/( 4*Delta lam* Delta lon ) sym = -1. call stencmd(zigm2_glb,cs,nmlon0,nmlat0,sym,cee,3) ! ! Insert RHS in finest stencil (formerly sub rths): do j = 1,nmlat0 jj = nmlath-nmlat0+j do i = 1,nmlon0 c0(i,j,10) = rhs_glb(i,jj) enddo ! i = 1,nmlon0 enddo ! j = 1,nmlat0 c0(nmlonp1,1,10) = c0(1,1,10) ! ! Set boundary condition at the pole: call edges(c0,nmlon0,nmlat0) call edges(c1,nmlon1,nmlat1) call edges(c2,nmlon2,nmlat2) call edges(c3,nmlon3,nmlat3) call edges(c4,nmlon4,nmlat4) call edges(cofum,nmlon0,nmlat0) ! ! Divide stencils by cos(lam_0) (not rhs): call divide(c0,nmlon0,nmlat0,nmlon0,nmlat0,cs,1) call divide(c1,nmlon1,nmlat1,nmlon0,nmlat0,cs,1) call divide(c2,nmlon2,nmlat2,nmlon0,nmlat0,cs,1) call divide(c3,nmlon3,nmlat3,nmlon0,nmlat0,cs,1) call divide(c4,nmlon4,nmlat4,nmlon0,nmlat0,cs,1) call divide(cofum,nmlon0,nmlat0,nmlon0,nmlat0,cs,0) ! ! Set value of solution to 1. at pole: do i=1,nmlon0 c0(i,nmlat0,10) = 1. enddo ! #if defined(INTERCOMM) || defined(CISMAH) ! ! 3/20/14 btf: cism coupling not yet ready in pdynamo (however, ! this part will not have to be parallelized, since we are in ! global domain here anyway). ! ! Convert LFM geographic potential to geomagnetic coordinates ! call cism_pot2mag c do jj=1,nmlat c do i=1,nmlonp1 c write(6,*)phihm(i,jj),i,jj c enddo c enddo #endif ! ! Modify stencils and RHS so that the NH high lat potential is inserted at ! high latitude. The SH high lat potential will be added back later. ! pfrac = fraction of dynamo in solution in the NH. = 1 low lat, = 0 hi lat ! cons_module: crit(1)=15, crit(2)=30 deg colats, or hi-lat > 75 deg, ! dynamo < 60 deg, and combination between 60-75 mag lat. ! The dynamo is symmetric about the magnetic equator, but the high latitude ! is anti-symmetric in both hemispheres. However, since Mudpack uses the ! NH potential pattern, then the SH potential pattern must be added ! back into the 2-D phim before the call threed, and before it is ! transformed to geographic coordinates. ! ! (In original dynamo, this was conditional on .not.mod_heelis) ! ncc = 1 nmaglon = nmlon0 nmaglat = nmlat0 do n=1,5 call stenmd(nmaglon,nmaglat,cee(ncc),phihm(1,nmlat0),pfrac) ncc = ncc+9*nmaglon*nmaglat if (n==1) ncc = ncc+nmaglon*nmaglat ! rhs is in 10th slot nmaglon = (nmaglon+1)/2 nmaglat = (nmaglat+1)/2 enddo ! n=1,5 ! ! (note nmlat0==nmlath, and nmlon0==nmlon+1) ! real :: c0_plt(nmlonp1,nmlat,10) ! real :: c0(nmlon0,nmlat0,10), ! c0_plt = 0. do n=1,10 c0_plt(:,1:nmlath,n) = c0(:,:,n) ! do j=1,nmlath ! write(6,"('n=',i2,' j=',i3,' c0(1,j,n)=',e15.7, ! | ' c0(nmlon:nmlonp1,j,n)=',2e15.7)") ! | n,j,c0_plt(1,j,n),c0_plt(nmlon:nmlonp1,j,n) ! enddo enddo ! call addfld('C0_01',' ',' ',c0_plt(:,:,1), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_02',' ',' ',c0_plt(:,:,2), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_03',' ',' ',c0_plt(:,:,3), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_04',' ',' ',c0_plt(:,:,4), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_05',' ',' ',c0_plt(:,:,5), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_06',' ',' ',c0_plt(:,:,6), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_07',' ',' ',c0_plt(:,:,7), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_08',' ',' ',c0_plt(:,:,8), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_09',' ',' ',c0_plt(:,:,9), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('C0_10',' ',' ',c0_plt(:,:,10), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) cofum_plt = 0. do n=1,9 cofum_plt(:,1:nmlath,n) = cofum(:,:,n) enddo ! call addfld('COFUM_01',' ',' ',cofum_plt(:,:,1), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('COFUM_02',' ',' ',cofum_plt(:,:,2), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('COFUM_03',' ',' ',cofum_plt(:,:,3), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('COFUM_04',' ',' ',cofum_plt(:,:,4), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('COFUM_05',' ',' ',cofum_plt(:,:,5), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('COFUM_06',' ',' ',cofum_plt(:,:,6), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('COFUM_07',' ',' ',cofum_plt(:,:,7), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('COFUM_08',' ',' ',cofum_plt(:,:,8), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) ! call addfld('COFUM_09',' ',' ',cofum_plt(:,:,9), ! | 'mlon',1,nmlonp1,'mlat',1,nmlat,0) end subroutine stencils !----------------------------------------------------------------------- subroutine clearcee(cee,nlon0,nlat0) ! ! Zero C arrays for stencil coefficients. ! Cee will contain: ! c0(nmlon0,nmlat0,10), c1(nmlon1,nmlat1,9), c2(nmlon2,nmlat2,9), ! c3(nmlon3,nmlat3,9), c4(nmlon4,nmlat4,9) ! ! Args: integer,intent(in) :: nlon0,nlat0 real,intent(out) :: cee(*) ! ! Local: integer :: nlon,nlat,n,m,i real :: time0,time1 ! if (dynamo_timing) call timer(time0,time1,"clearcee",0,0) ! ! Compute total size of cee nlon = nlon0 nlat = nlat0 n = 0 do m=1,5 ! 5 resolution levels n = n+nlon*nlat nlon = (nlon+1)/2 nlat = (nlat+1)/2 enddo ! m=1,5 ! 5 resolution levels n = 9*n+nlon0*nlat0 ! ! Clear cee: do i=1,n cee(i) = 0. enddo end subroutine clearcee !----------------------------------------------------------------------- subroutine stencmd(zigm,cs,nlon0,nlat0,sym,cee,ncoef) ! ! Calculate contribution fo 3 by 3 stencil from coefficient zigm ! at each grid point and level. ! ! Args: integer,intent(in) :: | nlon0, ! longitude dimension of finest grid level | nlat0, ! latitude dimension of finest grid level | ncoef ! integer identifier of coefficient real,intent(in) :: | zigm(nlon0,nlat0), ! coefficients (nlon0+1/2,(nlat0+1)/2) | sym, ! 1. if zigm symmetric w.r.t. equator ! -1. for antisymmetry | cs(nlat0) real,intent(inout) :: | cee(*) ! output stencil array consisting of c0,c1,c2,c3,c4 ! ! Local: integer :: nc,nlon,nlat,n real :: wkarray(-15:nmlon0+16,nmlat0) ! ! Perform half-way interpolation and extend zigm in wkarray: ! call htrpex(zigm,nlon0,nlat0,sym,wkarray) ! ! Calculate contribution to stencil for each grid point and level: ! nc = 1 nlon = nlon0 nlat = nlat0 ! ! Calculate modified and unmodified stencil on finest grid ! call cnmmod(nlon0,nlat0,nlon,nlat,cee(nc),ncoef,wkarray,cofum) ! if (dynamo_timing) call timer(time0,time1,"stencmd",0,0) ! ! Stencils on other grid levels remain the same. nc = nc+10*nlon*nlat nlon = (nlon+1)/2 nlat = (nlat+1)/2 ! do n=2,5 call cnm(nlon0,nlat0,nlon,nlat,cee(nc),ncoef,wkarray) nc = nc+9*nlon*nlat if (n==1) nc = nc+nlon*nlat nlon = (nlon+1)/2 nlat = (nlat+1)/2 enddo ! n=1,5 end subroutine stencmd !----------------------------------------------------------------------- subroutine htrpex(coeff,nmlon0,nmlat0,sym,wkarray) ! ! Perform half-way interpolation on array coeff and extend over 16 grid ! points. Result returned in wkarray. ! ! Args: integer,intent(in) :: nmlon0,nmlat0 real,intent(in) :: coeff(nmlon0,nmlat0),sym real,intent(out) :: wkarray(-15:nmlon0+16,nmlat0) ! ! Local: integer :: i,j,jj ! ! Copy coeff into positions in wkarray: do j=1,nmlat0 jj = nmlat0-j+1 do i=1,nmlon0 wkarray(i,j) = sym*coeff(i,jj) enddo ! i=1,nmlon0 enddo ! j=1,nmlat0 ! ! Extend over 32 grid spaces to allow for a total of 5 grid levels: do i=1,16 do j=1,nmlat0 wkarray(1-i,j) = wkarray(nmlon0-i,j) wkarray(nmlon0+i,j) = wkarray(1+i,j) enddo ! j=1,nmlat0 enddo ! i=1,16 end subroutine htrpex !----------------------------------------------------------------------- subroutine cnm(nlon0,nlat0,nlon,nlat,c,ncoef,wkarray) ! ! Compute contribution to stencil from zigm(ncoef) on grid nlon by nlat, ! Finest grid is nlon0 by nlat0. ! ! Args: integer,intent(in) :: | nlon0,nlat0, ! finest grid dimensions | nlon,nlat ! output grid dimensions real,intent(in) :: wkarray(-15:nmlon0+16,nmlat0) ! ! ncoef: integer id of coefficient: ! ncoef = 1 for zigm11 ! ncoef = 2 for zigm12 (=zigmc+zigm2) ! ncoef = 3 for zigm21 (=zigmc-zigm2) ! ncoef = 4 for zigm22 ! integer,intent(in) :: ncoef real,intent(inout) :: | c(nlon,nlat,*) ! output array for grid point stencils at ! resolution nlon x nlat ! ! Local: integer :: i,j,nint,i0,j0 real,parameter :: pi=3.141592654 real :: wk(nlon0,3) real :: time0,time1 ! if (dynamo_timing) call timer(time0,time1,"cnm",0,0) ! ! Compute separation of grid points of resolution nlon x nlat within ! grid of resolution nlon0,nlat0. Evaluate dlon and dlat, grid spacing ! of nlon x nlat. ! nint = (nlon0-1)/(nlon-1) ! ! Scan wkarray nlon x nlat calculating and adding contributions to stencil ! from zigm(ncoef) i0 = 1-nint j0 = 1-nint ! ! zigm11: ! am 2001-6-27 include boundary condition at equator if (ncoef==1) then do j = 1,nlat-1 do i = 1,nlon c(i,j,1) = c(i,j,1)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i+1)*nint,j0+j*nint)) c(i,j,5) = c(i,j,5)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) c(i,j,9) = c(i,j,9)-.5*(wkarray(i0+(i+1)*nint,j0+j*nint)+ | 2.*wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) enddo ! i = 1,nlon enddo ! j = 2,nlat-1 ! ! zigm12 (=zigmc+zigm2) elseif (ncoef==2) then do j = 2,nlat-1 do i = 1,nlon c(i,j,2) = c(i,j,2)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i+1)*nint,j0+j*nint)) c(i,j,4) = c(i,j,4)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) c(i,j,6) = c(i,j,6)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) c(i,j,8) = c(i,j,8)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i+1)*nint,j0+j*nint)) wk(i,1) = .5*(wkarray(i0+(i+1)*nint,j0+j*nint)- | wkarray(i0+(i-1)*nint,j0+j*nint)) wk(i,2) = (c(i,j,3)+wk(i,1))*(c(i,j,7)-wk(i,1)) wk(i,3) = sign(wk(i,1),c(i,j,3)+c(i,j,7)) if (wk(i,2) >= 0.) wk(i,3) = 0. c(i,j,3) = c(i,j,3)+wk(i,1)+wk(i,3) c(i,j,7) = c(i,j,7)-wk(i,1)+wk(i,3) c(i,j,9) = c(i,j,9)-2.*wk(i,3) enddo ! i = 1,nlon enddo ! j = 2,nlat-1 ! ! zigm21 (=zigmc-zigm2) elseif (ncoef==3) then do j = 2,nlat-1 do i = 1,nlon c(i,j,2) = c(i,j,2)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,4) = c(i,j,4)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,6) = c(i,j,6)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j-1)*nint)) c(i,j,8) = c(i,j,8)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j-1)*nint)) wk(i,1) = .5*(wkarray(i0+i*nint,j0+(j+1)*nint)- | wkarray(i0+i*nint,j0+(j-1)*nint)) wk(i,2) = (c(i,j,1)+wk(i,1))*(c(i,j,5)-wk(i,1)) wk(i,3) = sign(wk(i,1),c(i,j,1)+c(i,j,5)) if (wk(i,2) >= 0.) wk(i,3) = 0. c(i,j,1) = c(i,j,1)+wk(i,1)+wk(i,3) c(i,j,5) = c(i,j,5)-wk(i,1)+wk(i,3) c(i,j,9) = c(i,j,9)-2.*wk(i,3) enddo ! i = 1,nlon enddo ! j = 2,nlat-1 ! ! Low latitude boundary condition: j = 1 do i=1,nlon c(i,j,2) = c(i,j,2)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,4) = c(i,j,4)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) wk(i,1) = .5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) wk(i,2) = (c(i,j,1)+wk(i,1))*(c(i,j,5)-wk(i,1)) wk(i,3) = sign(wk(i,1),c(i,j,1)+c(i,j,5)) if (wk(i,2) >= 0.) wk(i,3) = 0. c(i,j,1) = c(i,j,1)+wk(i,1)+wk(i,3) c(i,j,5) = c(i,j,5)-wk(i,1)+wk(i,3) c(i,j,9) = c(i,j,9)-2.*wk(i,3) enddo ! i=1,nlon ! ! zigm22: elseif (ncoef==4) then do j = 2,nlat-1 do i = 1,nlon c(i,j,3) = c(i,j,3)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,7) = c(i,j,7)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j-1)*nint)) c(i,j,9) = c(i,j,9)-.5*(wkarray(i0+i*nint,j0+(j-1)*nint) | +2.*wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j+1)*nint)) enddo ! i = 1,nlon enddo ! j = 2,nlat-1 ! ! Low latitude boundary condition: j = 1 do i=1,nlon c(i,j,3) = c(i,j,3)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,9) = c(i,j,9)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) enddo ! i=1,nlon endif ! ncoef end subroutine cnm !----------------------------------------------------------------------- subroutine cnmmod(nlon0,nlat0,nlon,nlat,c,ncoef,wkarray,cofum) ! ! Compute contribution to stencil from zigm(ncoef) on grid nlon by nlat, ! Finest grid is nlon0 by nlat0. ! ! Args: integer,intent(in) :: | nlon0,nlat0, ! finest grid dimensions | nlon,nlat ! output grid dimensions real,intent(in) :: wkarray(-15:nmlon0+16,nmlat0) real,dimension(nmlon0,nmlat0,9),intent(inout) :: cofum ! ! ncoef: integer id of coefficient: ! ncoef = 1 for zigm11 ! ncoef = 2 for zigm12 (=zigmc+zigm2) ! ncoef = 3 for zigm21 (=zigmc-zigm2) ! ncoef = 4 for zigm22 ! integer,intent(in) :: ncoef real,intent(inout) :: | c(nlon,nlat,*) ! output array for grid point stencils at ! resolution nlon x nlat ! ! Local: integer :: i,j,nint,i0,j0 real,parameter :: pi=3.141592654 real :: wk(nlon0,3) ! ! Compute separation of grid points of resolution nlon x nlat within ! grid of resolution nlon0,nlat0. Evaluate dlon and dlat, grid spacing ! of nlon x nlat. ! nint = (nlon0-1)/(nlon-1) ! ! Scan wkarray nlon x nlat calculating and adding contributions to stencil ! from zigm(ncoef) i0 = 1-nint j0 = 1-nint ! ! zigm11: ! am 2001-6-27 include boundary condition at equator if (ncoef==1) then do j = 1,nlat-1 do i = 1,nlon c(i,j,1) = c(i,j,1)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i+1)*nint,j0+j*nint)) c(i,j,5) = c(i,j,5)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) c(i,j,9) = c(i,j,9)-.5*(wkarray(i0+(i+1)*nint,j0+j*nint)+ | 2.*wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) ! ! Unmodified: cofum(i,j,1) = cofum(i,j,1)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+(i+1)*nint,j0+j*nint)) cofum(i,j,5) = cofum(i,j,5)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+(i-1)*nint,j0+j*nint)) cofum(i,j,9) = cofum(i,j,9)- | .5*(wkarray(i0+(i+1)*nint,j0+j*nint)+ | 2.*wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) enddo ! i = 1,nlon enddo ! j = 2,nlat-1 ! ! zigm12 (=zigmc+zigm2) elseif (ncoef==2) then do j = 2,nlat-1 do i = 1,nlon c(i,j,2) = c(i,j,2)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i+1)*nint,j0+j*nint)) c(i,j,4) = c(i,j,4)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) c(i,j,6) = c(i,j,6)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i-1)*nint,j0+j*nint)) c(i,j,8) = c(i,j,8)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+(i+1)*nint,j0+j*nint)) wk(i,1) = .5*(wkarray(i0+(i+1)*nint,j0+j*nint)- | wkarray(i0+(i-1)*nint,j0+j*nint)) ! ! Unmodified: cofum(i,j,2) = c(i,j,2) cofum(i,j,4) = c(i,j,4) cofum(i,j,6) = c(i,j,6) cofum(i,j,8) = c(i,j,8) cofum(i,j,3) = cofum(i,j,3)+wk(i,1) cofum(i,j,7) = cofum(i,j,7)-wk(i,1) ! wk(i,2) = (c(i,j,3)+wk(i,1))*(c(i,j,7)-wk(i,1)) wk(i,3) = sign(wk(i,1),c(i,j,3)+c(i,j,7)) if (wk(i,2) >= 0.) wk(i,3) = 0. c(i,j,3) = c(i,j,3)+wk(i,1)+wk(i,3) c(i,j,7) = c(i,j,7)-wk(i,1)+wk(i,3) c(i,j,9) = c(i,j,9)-2.*wk(i,3) enddo ! i = 1,nlon enddo ! j = 2,nlat-1 ! ! zigm21 (=zigmc-zigm2) elseif (ncoef==3) then do j = 2,nlat-1 do i = 1,nlon c(i,j,2) = c(i,j,2)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,4) = c(i,j,4)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,6) = c(i,j,6)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j-1)*nint)) c(i,j,8) = c(i,j,8)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j-1)*nint)) wk(i,1) = .5*(wkarray(i0+i*nint,j0+(j+1)*nint)- | wkarray(i0+i*nint,j0+(j-1)*nint)) ! ! Unmodified: cofum(i,j,2) = c(i,j,2) cofum(i,j,4) = c(i,j,4) cofum(i,j,6) = c(i,j,6) cofum(i,j,8) = c(i,j,8) cofum(i,j,1) = cofum(i,j,1)+wk(i,1) cofum(i,j,5) = cofum(i,j,5)-wk(i,1) ! wk(i,2) = (c(i,j,1)+wk(i,1))*(c(i,j,5)-wk(i,1)) wk(i,3) = sign(wk(i,1),c(i,j,1)+c(i,j,5)) if (wk(i,2) >= 0.) wk(i,3) = 0. c(i,j,1) = c(i,j,1)+wk(i,1)+wk(i,3) c(i,j,5) = c(i,j,5)-wk(i,1)+wk(i,3) c(i,j,9) = c(i,j,9)-2.*wk(i,3) enddo ! i = 1,nlon enddo ! j = 2,nlat-1 ! ! Low latitude boundary condition: j = 1 do i=1,nlon c(i,j,2) = c(i,j,2)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,4) = c(i,j,4)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) wk(i,1) = .5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) cofum(i,j,2) = c(i,j,2) cofum(i,j,4) = c(i,j,4) cofum(i,j,1) = cofum(i,j,1)+wk(i,1) cofum(i,j,5) = cofum(i,j,5)-wk(i,1) wk(i,2) = (c(i,j,1)+wk(i,1))*(c(i,j,5)-wk(i,1)) wk(i,3) = sign(wk(i,1),c(i,j,1)+c(i,j,5)) if (wk(i,2) >= 0.) wk(i,3) = 0. c(i,j,1) = c(i,j,1)+wk(i,1)+wk(i,3) c(i,j,5) = c(i,j,5)-wk(i,1)+wk(i,3) c(i,j,9) = c(i,j,9)-2.*wk(i,3) enddo ! i=1,nlon ! ! zigm22: elseif (ncoef==4) then do j = 2,nlat-1 do i = 1,nlon c(i,j,3) = c(i,j,3)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,7) = c(i,j,7)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j-1)*nint)) c(i,j,9) = c(i,j,9)-.5*(wkarray(i0+i*nint,j0+(j-1)*nint) | +2.*wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j+1)*nint)) ! ! Unmodified: cofum(i,j,3) = cofum(i,j,3)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j+1)*nint)) cofum(i,j,7) = cofum(i,j,7)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j-1)*nint)) cofum(i,j,9) = cofum(i,j,9)-.5*(wkarray(i0+i*nint,j0+(j-1)* | nint)+2.*wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) enddo ! i = 1,nlon enddo ! j = 2,nlat-1 ! ! Low latitude boundary condition: j = 1 do i=1,nlon c(i,j,3) = c(i,j,3)+.5*(wkarray(i0+i*nint,j0+j*nint) | +wkarray(i0+i*nint,j0+(j+1)*nint)) c(i,j,9) = c(i,j,9)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) cofum(i,j,3) = cofum(i,j,3)+.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) cofum(i,j,9) = cofum(i,j,9)-.5*(wkarray(i0+i*nint,j0+j*nint)+ | wkarray(i0+i*nint,j0+(j+1)*nint)) enddo ! i=1,nlon endif ! ncoef end subroutine cnmmod !----------------------------------------------------------------------- subroutine ceee(cee,nx,ny,cf) ! ! Called from mudpack solvers to transfer coefficients. ! save integer,intent(in) :: nx,ny real,intent(in) :: cee(nx,ny,*) real,intent(out) :: cf(nx,ny,*) integer :: i,j,n real :: time0,time1 do n = 1,9 do j = 1,ny do i = 1,nx cf(i,j,n) = cee(i,j,n) enddo ! write(6,"('ceee: n=',i2,' j=',i3,' nx=',i3,' ny=',i3, ! | ' cf(:,j,n)=',/,(6e12.4))") n,j,nx,ny,cf(:,j,n) enddo enddo end subroutine ceee !----------------------------------------------------------------------- subroutine edges(c,nlon,nlat) ! ! Insert equatorial and polar boundary conditions in stencil c(nlon,nlat,9) ! ! Args: integer,intent(in) :: nlon,nlat real,intent(out) :: c(nlon,nlat,*) ! ! Local: integer :: n,i do n=1,8 do i=1,nlon c(i,nlat,n) = 0. enddo enddo do i=1,nlon c(i,nlat,9) = 1. enddo end subroutine edges !----------------------------------------------------------------------- subroutine divide(c,nlon,nlat,nlon0,nlat0,cs,igrid) ! ! Divide stencil C by cos(theta(i,j)) ! ! Args: integer,intent(in) :: nlon,nlat,nlon0,nlat0,igrid real,intent(in) :: cs(*) real,intent(out) :: c(nlon,nlat,*) ! ! Local: integer :: nint,j0,n,j,i ! nint = (nlon0-1)/(nlon-1) j0 = 1-nint do n = 1,9 do j = 1,nlat-1 do i = 1,nlon c(i,j,n) = c(i,j,n)/(cs(j0+j*nint)*nint**2) enddo ! i = 1,nlon enddo ! j = 1,nlat-1 enddo ! n = 1,9 ! ! This is unnecessary since cs(1)==1. ! (this was in serial dynamo of timegcm, but not in tiegcm). ! ! if (nint==1.and.igrid > 0) then ! do i = 1,nlon ! c(i,1,10) = c(i,1,10)/cs(1) ! enddo ! i = 1,nlon ! endif end subroutine divide !----------------------------------------------------------------------- subroutine stenmd(inlon,inlat,c,phihm,pfrac) use cons_module,only: dlatm,dtr ! ! Modify stencil to set potential to heelis value within auroral circle. ! ! Args: integer,intent(in) :: inlon,inlat real,intent(inout) :: c(inlon,inlat,*) real,dimension(nmlon0,nmlat0),intent(in) :: | phihm, ! heelis potential (from subs potm, flwv32) | pfrac ! fractional presence of dynamo (from sub colath) ! ! Local: integer :: nint,i0,j0,i,j,n,jj ! ! Compute separation of grid points for this resolution: nint = (nmlon0-1)/(inlon-1) i0 = 1-nint j0 = 1-nint ! ! If nint==1, then we are at the highest resolution. ! Correct RHS, which is in c(10) ! if (nint==1) then do j=1,inlat do i=1,inlon c(i,j,10) = pfrac(i,j)*c(i,j,10)+(1.-pfrac(i,j))*c(i,j,9)* | (dlatm/(10.*dtr))**2*phihm(i,j) enddo ! i=1,inlon enddo ! j=1,inlat endif ! ! Modify stencil, c(i,j,n),n=1,9: ! if (nint==1) then do j=1,inlat jj = j0+j*nint do n = 1,8 do i = 1,inlon c(i,j,n) = c(i,j,n)*pfrac(i0+i*nint,jj) cofum(i,j,n) = cofum(i,j,n)*pfrac(i0+i*nint,jj) enddo ! i = 1,inlon enddo ! n = 1,8 do i = 1,inlon c(i,j,9) = c(i,j,9)*pfrac(i0+i*nint,jj)+ | (1.-pfrac(i0+i*nint,jj))*c(i,j,9)* | (dlatm*float(nint)/(10.*dtr))**2 cofum(i,j,9) =cofum(i,j,9)*pfrac(i0+i*nint,jj)+ | (1.-pfrac(i0+i*nint,jj))*cofum(i,j,9)* | (dlatm*float(nint)/(10.*dtr))**2 enddo ! i = 1,inlon enddo ! j=1,inlat else ! nint /= 1 do j=1,inlat jj = j0+j*nint do n = 1,8 do i = 1,inlon c(i,j,n) = c(i,j,n)*pfrac(i0+i*nint,jj) enddo ! i = 1,inlon enddo ! n = 1,8 do i = 1,inlon c(i,j,9) = c(i,j,9)*pfrac(i0+i*nint,jj)+ | (1.-pfrac(i0+i*nint,jj))*c(i,j,9)* | (dlatm*float(nint)/(10.*dtr))**2 enddo ! i = 1,inlon enddo ! j=1,inlat endif ! nint end subroutine stenmd !----------------------------------------------------------------------- subroutine solver(cofum,c0) ! ! Call mudpack to solve PDE. Solution is returned in rim: ! real,dimension(nmlonp1,nmlat,2) :: rim ! Mudpack solvers: ! isolve = 0 org. mud v5. (modified stencils neq direct solution) ! isolve = 1 muh hybrid solver (no convergence => only as direct solver) ! isolve = 2 modified mud (residual calculated with unmodified stencils ! same solution as with direct solver, if same ! coefficient matrix is used) ! Only isolve==2 is supported in pdynamo. ! ! Args: real,dimension(nmlon0,nmlat0,9),intent(in) :: cofum real,intent(in) :: c0(nmlon0,nmlat0,10) ! ! Locals: integer :: i,j,jntl,ier,isolve real :: l2norm real,dimension(nmlon0,nmlat0) :: ressolv real,dimension(nmlon0,nmlat0,9) :: cofum_solv ! Module data: ! real,dimension(nmlonp1,nmlat,2) :: rim_glb ! pde solver output ! call addfld('RIM1_PRESOLV',' ','A/m',rim_glb(:,:,1),'mlon',1, ! | nmlonp1,'mlat',1,nmlat,0) ! call addfld('RIM2_PRESOLV',' ','A/m',rim_glb(:,:,2),'mlon',1, ! | nmlonp1,'mlat',1,nmlat,0) jntl = 0 ier = 0 isolve = 2 call mudmod(rim_glb,phisolv,jntl,isolve,ier)! solver in mudmod.F if (ier < 0 ) then ! not converged write(6,*) 'muh: use direct solver' call muh(rim_glb,jntl) ! solver in mud.F endif ! call addfld('RIM1_SOLV',' ','A/m',rim_glb(:,:,1),'mlon',1,nmlonp1, ! | 'mlat',1,nmlat,0) ! call addfld('RIM2_SOLV',' ','A/m',rim_glb(:,:,2),'mlon',1,nmlonp1, ! | 'mlat',1,nmlat,0) l2norm=0. ressolv = 0.0 do j = 1,nmlat0 do i = 1,nmlon0-1 cofum_solv(i,j,:)= cofum(i,j,:) ! ! fields: phisolv(0:nmlonp1,0:nmlat+1) ! 2d solution/ electric potential ! ressolv(i,j) = ( + cofum_solv(i,j,1)*phisolv(i+1,j)+ + cofum_solv(i,j,2)*phisolv(i+1,j+1)+ + cofum_solv(i,j,3)*phisolv(i,j+1)+ + cofum_solv(i,j,4)*phisolv(i-1,j+1)+ + cofum_solv(i,j,5)*phisolv(i-1,j)+ + cofum_solv(i,j,6)*phisolv(i-1,j-1)+ + cofum_solv(i,j,7)*phisolv(i,j-1)+ + cofum_solv(i,j,8)*phisolv(i+1,j-1)+ + cofum_solv(i,j,9)*phisolv(i,j)) ressolv(i,j) = c0(i,j,10)-ressolv(i,j) l2norm = l2norm + ressolv(i,j)*ressolv(i,j) end do end do ! write(6,*) 'L2norm (global root task) = ',l2norm end subroutine solver !----------------------------------------------------------------------- #if defined(INTERCOMM) || defined(CISMAH) subroutine cism_pot2mag ! ! This subroutine converts high latitude potential in geographic coordinate ! obtained form the CMIT M-I couplier to geomagnetci coordinate to be used ! in dynamo calculations ! ! Uses ! use params_module,only: nlonp1 use cism_coupling_module, only:gpotm use magfield_module,only: ig,jg,wt ! ! Local: ! integer :: jj do jj=mlat0,mlat1 call geo2mag(phihm(1,jj),gpotm(3:nlonp1+2,:), | ig,jg,wt,nlonp1,nmlonp1,nmlon,nmlat,jj) ! Periodic point phihm(nmlonp1,jj) = phihm(1,jj) enddo end subroutine cism_pot2mag !----------------------------------------------------------------------- subroutine geo2mag(fmag,fgeo,long,latg,wght,nlonp1_geo,nlonp1_mag, | nlon_mag,nlat_mag,lat) ! ! Transform field fgeo on geographic grid to geomagnetic grid using ! indices long,latg and weights wght. Return field fmag on magnetic ! grid. ! ! Args: integer,intent(in) :: nlonp1_geo,nlonp1_mag,nlon_mag,nlat_mag,lat integer,dimension(nlonp1_mag,nlat_mag),intent(in) :: long,latg real,intent(in) :: fgeo(nlonp1_geo,*),wght(4,nlonp1_mag,nlat_mag) real,intent(out) :: fmag(nlonp1_mag,*) ! ! Local: integer :: i ! do i=mlon0,mlon1 fmag(i,1) = | fgeo(long(i,lat) ,latg(i,lat) )*wght(1,i,lat)+ | fgeo(long(i,lat)+1,latg(i,lat) )*wght(2,i,lat)+ | fgeo(long(i,lat)+1,latg(i,lat)+1)*wght(3,i,lat)+ | fgeo(long(i,lat) ,latg(i,lat)+1)*wght(4,i,lat) ! if (iprint > 0) write(6,"('geo2mag: i=',i3,' lat=',i3,' long=', ! | i3,' latg=',i3,' wght=',4e12.4,' fgeo=',e12.4,' fmag=', ! | e12.4)") i,lat,long(i,lat),latg(i,lat),wght(:,i,lat), ! | fgeo(long(i,lat),latg(i,lat)),fmag(i,1) enddo end subroutine geo2mag #endif !----------------------------------------------------------------------- end module pdynamo_module