module heelis use shr_kind_mod ,only: r8 => shr_kind_r8 ! 8-byte reals use maggrid ,only: nmlon,nmlonp1,nmlat,ylonm,ylatm use geogrid ,only: nlat,nlon use params ,only: dtr,rtd,pi,iulog use read_ncfile ,only: ctpoten use addfld_mod ,only: addfld ! routine to add field to output use namelist ,only: model_name ! ! phihm and pfrac are output of this module: ! use solve,only: & ! output nmlat0 ,& ! (nmlat+1)/2 phihm ,& ! high-latitude potential (nmlonp1,nmlat) pfrac ! NH fraction of potential (nmlonp1,nmlat0) implicit none save private ! ! Auroral parameters (taken from aurora.F of timegcm): ! (dimension of 2 is for south,north hemispheres) ! real(r8) :: & offa(2) ,& ! offset of oval towards 0 MLT relative to magnetic pole (rad) offc(2) ,& ! dskofa(2) ,& ! offset of oval in radians towards 18 MLT (f(By)) dskofc(2) ,& ! phin(2) ,& ! night convection entrance in MLT converted to radians (f(By)) phid(2) ,& ! dayside convection entrance in MLT converted to radians (f(By)) psim(2) ,& ! night convection entrance in MLT converted to radians (f(By)) psie(2) ,& ! pcen(2) ,& ! phidp0(2) ,& ! phidm0(2) ,& ! phinp0(2) ,& ! phinm0(2) ,& ! theta0(2) ,& ! convection reversal boundary in radians rr1(2) ! ! ! Critical angles (radians south,north) of transition from edynamo ! potential to Heelis potential (taken from timegcm cons.F) ! real(r8),parameter :: crit(2) = (/0.261799387_r8, 0.523598775_r8/) ! ! Sun's longitude in dipole coordinates (see sub sunloc) ! real(r8),allocatable :: sunlons(:) ! (nlat) ! public heelis_model,pfrac,phihm contains !----------------------------------------------------------------------- subroutine heelis_model(iday,secs) ! ! Driver for Heelis empirical model to calculate high-latitude potential. ! ! Args: integer,intent(in) :: iday ! day of year real(r8),intent(in) :: secs ! ut in seconds ! ! Local: real(r8),dimension(nmlonp1,nmlat) :: pfrac_plt ! ! Set auroral parameters: ! call heelis_init ! ! Get sun's location: ! call sunloc(iday,secs) ! ! Calculate pfrac: ! call colath ! ! Calculate the heelis potential phihm in geomagnetic coordinates: ! (potm calls sub flwv32) ! call potm ! pfrac_plt = 0._r8 ! pfrac_plt(:,1:nmlat0) = pfrac(:,:) ! call addfld('PHIHM',' ',' ',phihm,'mlon',1,nmlonp1,& ! 'mlat',1,nmlat,0) ! call addfld('PFRAC',' ',' ',pfrac_plt,'mlon',1,nmlonp1,& ! 'mlat',1,nmlat,0) end subroutine heelis_model !----------------------------------------------------------------------- subroutine heelis_init ! ! Auroral parameters for Heelis (taken from aurora.F of timegcm): ! ! This is called at every timestep because ctpoten may change with time. ! Time-dependent ctpoten (kV) is read from TIMEGCM and WACCM input files, ! unless it was provided as a constant by the user via namelist (namelist.F90). ! offa(:) = 1._r8*dtr offc(:) = 1._r8*dtr dskofa(:) = 0._r8 dskofc(:) = 0._r8 phin(:) = 180._r8*dtr phid(:) = 0._r8 psim(:) = 0.5_r8 * ctpoten * 1000._r8 psie(:) = -0.5_r8 * ctpoten * 1000._r8 pcen(:) = 0._r8 phidp0(:) = 90._r8*dtr phidm0(:) = 90._r8*dtr phinp0(:) = 90._r8*dtr phinm0(:) = 90._r8*dtr rr1(:) = -2.6_r8 theta0(:) = (-3.80_r8+8.48_r8*(ctpoten**0.1875_r8))*dtr ! write(iulog,"('heelis_init: ctpoten=',e12.4,' theta0=',e12.4)") ctpoten,theta0(1) end subroutine heelis_init !----------------------------------------------------------------------- subroutine colath ! ! Local: integer :: j,i real(r8),dimension(nmlonp1,nmlat0) :: colatc real(r8) :: sinlat,coslat,aslonc,ofdc,cosofc,sinofc ! ! offc(2), dskofc(2) are for northern hemisphere aurora ofdc = sqrt(offc(2)**2+dskofc(2)**2) cosofc = cos(ofdc) sinofc = sin(ofdc) aslonc = asin(dskofc(2)/ofdc) ! ! Define colatc with northern convection circle coordinates ! sunlons(nlat): sun's longitude in dipole coordinates (see sub sunloc) ! do j=1,nmlat0 sinlat = sin(abs(ylatm(j+nmlat0-1))) coslat = cos( ylatm(j+nmlat0-1)) do i=1,nmlonp1 colatc(i,j) = cos(ylonm(i)-sunlons(1)+aslonc) colatc(i,j) = acos(cosofc*sinlat-sinofc*coslat*colatc(i,j)) enddo ! i=1,nmlonp1 ! write(iulog,"('colath: j=',i3,' colatc(:,j)=',/,(6e12.4))") & ! j,colatc(:,j)*rtd ! ! Calculate fractional presence of dynamo equation at each northern ! hemisphere geomagnetic grid point. Output in pfrac(nmlonp1,nmlat0) ! do i=1,nmlonp1 pfrac(i,j) = (colatc(i,j)-crit(1))/(crit(2)-crit(1)) if (pfrac(i,j) < 0._r8) pfrac(i,j) = 0._r8 if (pfrac(i,j) >= 1._r8) pfrac(i,j) = 1._r8 enddo ! i=1,nmlonp1 ! write(iulog,"('colath: j=',i3,' pfrac(:,j)=',/,(6e12.4))") j,pfrac(:,j) enddo ! j=1,nmlat0 end subroutine colath !----------------------------------------------------------------------- subroutine potm use params ,only: pi_dyn ! pi used in dynamo calculations use maggrid,only: ylonm,ylatm ! magnetic grid lons, lats ! ! Calculate heelis potential in geomagnetic coordinates. ! ! Local: integer :: i,j real(r8),dimension(nmlon) :: dlat,dlon,ratio integer,dimension(nmlon) :: iflag ! ratio(:) = 1._r8 do j=1,nmlat iflag(:) = 1 ! must be updated at each j dlat(:) = ylatm(j) dlon(:) = ylonm(1:nmlon)-sunlons(1) ! ! flwv32 returns single-level Heelis potential in geomag coords: ! if (abs(ylatm(j)) > pi_dyn/6._r8) then call flwv32(dlat,dlon,ratio,iflag,nmlon,phihm(:,j),j) else phihm(1:nmlon,j) = 0._r8 endif enddo ! j=1,nmlat ! ! Periodic point: do j=1,nmlat phihm(nmlonp1,j) = phihm(1,j) enddo ! j=1,nmlat end subroutine potm !----------------------------------------------------------------------- subroutine flwv32(dlat,dlon,ratio,iflag,nmlon,poten,mlat) ! ! Calculate heelis potential at current magnetic latitude mlat. ! use params ,only: pi_dyn ! ! Args: integer,intent(in) :: mlat,nmlon integer,intent(inout) :: iflag(nmlon) real(r8),dimension(nmlon),intent(in) :: dlat,dlon,ratio real(r8),dimension(nmlon+1),intent(out) :: poten ! ! Local: integer :: i,n,ihem real(r8),parameter :: eps=1.e-6_r8 real(r8) :: & pi2,pih,sinthr1,psi(8),phirc,sinth0, & ofda,cosofa(2),sinofa(2),aslona(2), & ofdc,cosofc(2),sinofc(2),aslonc(2), & phdpmx(2),phnpmx(2),phnmmx(2),phdmmx(2) real(r8),dimension(nmlon) :: sinlat,coslat,sinlon,coslon,alon, & colat,wk1,wk2,wk3,phifun,phifn2 integer :: ifn(nmlon) real(r8) :: phi(nmlon,8) ! pi2 = 2.0_r8*pi_dyn pih = 0.5_r8*pi_dyn do n=1,2 ofda = sqrt(offa(n)**2+dskofa(n)**2) cosofa(n) = cos(ofda) sinofa(n) = sin(ofda) aslona(n) = asin(dskofa(n)/ofda) ! ofdc = sqrt(offc(n)**2+dskofc(n)**2) cosofc(n) = cos(ofdc) sinofc(n) = sin(ofdc) aslonc(n) = asin(dskofc(n)/ofdc) ! if (phin(n) < phid(n)) phin(n) = phin(n)+pi2 ! modifies aurora phin phdpmx(n) = .5_r8*min(pi,(phin(n)-phid(n))) phnpmx(n) = .5_r8*min(pi,(phid(n)-phin(n)+pi2)) phnmmx(n) = phdpmx(n) phdmmx(n) = phnpmx(n) enddo ! n=1,2 ! ! Set ihem=1,2 for South,North hemisphere: ! ihem = int(dlat(max0(1,nlon/2))*2._r8/3.1416_r8+2._r8) sinth0 = sin(theta0(ihem)) ! ! Average amie results show r1=-2.6 for 11.3 degrees ! (0.1972 rad) beyond theta0. ! sinthr1 = sin(theta0(ihem)+0.1972_r8) psi(1) = psie(ihem) psi(3) = psim(ihem) do n=2,4,2 psi(n) = psi(n-1) enddo ! n=2,4,2 do n=1,4 psi(n+4) = psi(n) enddo ! n=1,4 ! ! Transform to auroral circle coordinates: ! do i=1,nmlon sinlat(i) = sin(abs(dlat(i))) coslat(i) = cos(dlat(i)) sinlon(i) = sin(dlon(i)+aslonc(ihem)) coslon(i) = cos(dlon(i)+aslonc(ihem)) colat(i) = cosofc(ihem)*sinlat(i)-sinofc(ihem)*coslat(i)* & coslon(i) alon(i) = mod(atan2(sinlon(i)*coslat(i),sinlat(i)* & sinofc(ihem)+cosofc(ihem)*coslat(i)*coslon(i))- & aslonc(ihem)+3._r8*pi_dyn,pi2)-pi_dyn colat(i) = acos(colat(i))*sqrt(ratio(i)) ! ! Boundaries for longitudinal function: ! wk1(i) = ((colat(i)-theta0(ihem))/theta0(ihem))**2 phi(i,4)=phid(ihem)+eps-min(phidm0(ihem)+wk1(i)* & (pih-phidm0(ihem)),phdmmx(ihem)) phi(i,5)=phid(ihem)-eps+min(phidp0(ihem)+wk1(I)* & (pih-phidp0(ihem)),phdpmx(ihem)) phi(i,6)=phin(ihem)+eps-min(phinm0(ihem)+wk1(i)* & (pih-phinm0(ihem)),phnmmx(ihem)) phi(i,7)=phin(ihem)-eps+min(phinp0(ihem)+wk1(i)* & (pih-phinp0(ihem)),phnpmx(ihem)) phi(i,1)=phi(i,5)-pi2 phi(i,2)=phi(i,6)-pi2 phi(i,3)=phi(i,7)-pi2 phi(i,8)=phi(i,4)+pi2 phifun(i)=0._r8 phifn2(i) = 0._r8 if (colat(i)-theta0(ihem) >= 0._r8) then ifn(i) = 3 else ifn(i) = 2 endif if (iflag(i) == 1) iflag(i) = ifn(i) ! ! Add ring current rotation to potential (phirc) ! phirc = 0._r8 wk2(i) = mod(alon(i)+phirc+2._r8*pi2+pi_dyn,pi2)-pi_dyn wk3(i) = mod(alon(i)+phirc+3._r8*pi2,pi2)-pi_dyn enddo ! i=1,nmlon ! ! Longitudinal variation: ! do n=1,7 do i=1,nmlon phifun(i)=phifun(i)+.25_r8*(psi(n)+psi(n+1)+(psi(n)- & psi(n+1))*cos(mod(pi_dyn*(wk2(i)-phi(i,n))/(phi(i,n+1)- & phi(i,n)),pi2)))*(1._r8-sign(1._r8,(wk2(i)-phi(i,n))* & (wk2(i)-phi(i,n+1)))) phifn2(i)=phifn2(i)+.25_r8*(psi(n)+psi(n+1)+(psi(n)- & psi(n+1))*cos(mod(pi_dyn*(wk3(i)-phi(i,n))/(phi(i,n+1)- & phi(i,n)),pi2)))*(1._r8-sign(1._r8,(wk3(i)-phi(i,n))* & (wk3(i)-phi(i,n+1)))) enddo enddo ! ! Evaluate total potential: ! do i=1,nmlon if (iflag(i)==2) then poten(i) = (2._r8*(pcen(ihem)-phifun(i))+(phifun(i)-phifn2(i))* & 0.75_r8)*(colat(i)/theta0(ihem))**3 + & (1.5_r8*(phifun(i)+phifn2(i))-3._r8*pcen(ihem))*(colat(i)/ & theta0(ihem))**2 + 0.75_r8*(phifun(i)-phifn2(i))*(colat(i)/ & theta0(ihem)) + pcen(ihem) else poten(i) = phifun(i)*(max(sin(colat(i)),sinth0)/sinth0)**rr1(ihem)* & exp(7._r8*(1._r8-max(sin(colat(i)),sinthr1)/sinthr1)) endif enddo ! write(iulog,"(/'flwv32: mlat=',i2,' ihem=',i2)") mlat,ihem ! write(iulog,"(' theta0(ihem)=',e12.4,' pcen(ihem)=',e12.4,' rr1(ihem)=',e12.4)") & ! theta0(ihem),pcen(ihem),rr1(ihem) ! write(iulog,"(' sinth0=',e12.4,' sinthr1=',e12.4)") sinth0,sinthr1 ! write(iulog,"(' iflag=',/,(20i3))") iflag ! write(iulog,"(' dlat =',/,(6e12.4))") dlat ! write(iulog,"(' dlon =',/,(6e12.4))") dlon ! write(iulog,"(' colat=',/,(6e12.4))") colat ! write(iulog,"(' phifun=',/,(6e12.4))") phifun ! write(iulog,"(' phifn2=',/,(6e12.4))") phifn2 ! write(iulog,"(' poten =',/,(6e12.4))") poten end subroutine flwv32 !----------------------------------------------------------------------- subroutine sunloc(iday,secs) use apex,only: alonm ! real(r8),dimension(nlonp1,0:nlatp1) ! ! Given day of year and ut, return sun's longitudes in dipole coordinates ! in sunlons(nlat) ! ! Args: integer,intent(in) :: iday ! day of year real(r8),intent(in) :: secs ! ut in seconds ! ! Local: integer :: j,i,ii,isun,jsun real(r8) :: glats,glons,pisun,pjsun,sndlons,csdlons real(r8) :: dphi,dlamda real(r8),allocatable,save :: rlonm(:,:) ! (nlon+4,nlat) real(r8) :: r8_nlat, r8_nlon real(r8) :: r8_isun, r8_jsun logical,save :: first = .true. ! ! Sun's geographic coordinates: r8_nlat = real(nlat) r8_nlon = real(nlon) glats = asin(.398749_r8*sin(2._r8*pi*(iday-80)/365._r8)) glons = pi*(1._r8-2._r8*secs/86400._r8) dphi = pi/r8_nlat dlamda = 2._r8*pi/r8_nlon if (first) then allocate(rlonm(nlon+4,nlat)) do j=1,nlat do i=1,nlon ii = i+2 rlonm(ii,j) = alonm(i,j) enddo enddo do i=1,2 do j=1,nlat rlonm(i,j) = rlonm(i+nlon,j) rlonm (i+nlon+2,j) = rlonm (i+2,j) enddo enddo allocate(sunlons(nlat)) endif ! first call pisun = (glons+pi)/dlamda+1._r8 pjsun = (glats+.5_r8*(pi-dphi))/dphi+1._r8 isun = int(pisun) jsun = int(pjsun) r8_isun = real(isun) r8_jsun = real(jsun) pisun = pisun-r8_isun pjsun = pjsun-r8_jsun sndlons = (1._r8-pisun)*(1._r8-pjsun)*sin(rlonm(isun+2,jsun))+ & pisun*(1._r8-pjsun) *sin(rlonm(isun+3,jsun))+ & pisun*pjsun *sin(rlonm(isun+3,jsun+1))+ & (1._r8-pisun)*pjsun *sin(rlonm(isun+2,jsun+1)) csdlons = (1._r8-pisun)*(1._r8-pjsun)*cos(rlonm(isun+2,jsun))+ & pisun*(1._r8-pjsun) *cos(rlonm(isun+3,jsun))+ & pisun*pjsun *cos(rlonm(isun+3,jsun+1))+ & (1._r8-pisun)*pjsun *cos(rlonm(isun+2,jsun+1)) sunlons(1) = atan2(sndlons,csdlons) do j = 2,nlat sunlons(j) = sunlons(1) enddo if (first) first = .false. end subroutine sunloc !----------------------------------------------------------------------- end module heelis