61,62c61 < | pi, ! 4.*atan(1.) < | crit ! may be changed from default if amie run --- > | pi ! 4.*atan(1.) 65,67d63 < use amie_module,only: crad,phida,ekvg,efxg,hpi_sh_amie, < | hpi_nh_amie,pcp_sh_amie,pcp_nh_amie < use input_module,only: iamie 154d149 < | rradp(2), ! radius of polar-cap circle in radians for drizls 214,215d208 < use hist_module,only: modeltime < use init_module,only: iyear 223a217 > integer :: ifcalccloc 229c223 < real :: arad(2),hp_amie,c25,c35 --- > real :: arad(2) 288,302c282,290 < ! Recalculate power and ctpoten if an amie run: < if (iamie > 0) then < ! hp_amie = max(hpi_sh_amie,hpi_nh_amie) < hp_amie = 0.5*(hpi_sh_amie+hpi_nh_amie) < write (6,"(1x,'old,amie power S,N =',4e12.4)") power,hp_amie, < | hpi_sh_amie,hpi_nh_amie < power = hp_amie < ! Assign ctpoten as average of cps,n in kV of Polar Cap Potential drop < write (6,"(1x,'old,amie CP SH,NH =',3e12.4)") ctpoten, < | pcp_sh_amie,pcp_nh_amie < ctpoten = 0.5 * (pcp_sh_amie + pcp_nh_amie) < write (6,"(1x,'theta0,offc,dskofc,phid,phin =',10e12.4)") < | theta0/dtr,offc/dtr,dskofc/dtr,phid/(dtr*15)-12, < | phin/(dtr*15)-12 < endif --- > ! If AMIE, CMIT, Weimer or other provide power and/or ctpoten, they should be called > ! before this, or should reset values when they are called > > ! Set flag to use offc, dskofc, theta0, etc from calccloc > ifcalccloc = 0 > ! CMIT provides power and ctpoten, and calccloc provides the offsets, radii, and phid,n > #if defined(INTERCOMM) || defined(CISMAH) > !! ifcalccloc = 1 > #endif 308,310d295 < < rhp = 14.20 + 0.96*plevel < rcp = -0.43 + 9.69 * (ctpoten**0.1875) 337,371c322,323 < ! 08/12 bae: Already have these parameters from calccloc if an amie run: < if (iamie .eq. 0) then < theta0(isouth) = (-3.80+8.48*(ctpoten**0.1875))*dtr < theta0(inorth) = theta0(isouth) < offa(isouth) = 1.0*dtr < offa(inorth) = 1.0*dtr < dskofa(isouth) = 0. < dskofa(inorth) = 0. < ! < arad(isouth) = max(rcp,rhp) < arad(inorth) = max(rcp,rhp) < rrad(isouth) = arad(isouth)*dtr < rrad(inorth) = arad(inorth)*dtr < < else < < ! 08/12 bae: Define rradp for drizzle if amie run < ! rradp(isouth) = crad(isouth)-1.5*dtr < ! rradp(inorth) = crad(inorth)-2.5*dtr < rradp(isouth) = crad(isouth) + h2*dtr < rradp(inorth) = crad(inorth) + h2*dtr < if (rradp(isouth) > 29.*dtr) rradp(isouth) = 29.*dtr < if (rradp(inorth) > 29.*dtr) rradp(inorth) = 29.*dtr < ! if (rradp(isouth) > 31.*dtr) rradp(isouth) = 31.*dtr < ! if (rradp(inorth) > 31.*dtr) rradp(inorth) = 31.*dtr < ! Redefine rrad (araddeg = craddeg + 2 in calccloc) < rrad(isouth) = rradp(isouth) + 0.5*dtr < rrad(inorth) = rradp(inorth) + 0.5*dtr < arad(isouth) = rrad(isouth)/dtr < arad(inorth) = rrad(inorth)/dtr < ! 08/12 bae: Revise crit to come from theta0 in colath < call colath < write(6,"('crad =',f6.2,' rrad=',f6.2,' rradp=',f6.2, < | ' pcp_nh_amie=',f8.1,' crit(2)=',2f6.2)") crad(1)/dtr, < | rrad(1)/dtr,rradp(1)/dtr,pcp_nh_amie,crit/dtr --- > rhp = 14.20 + 0.96*plevel > rcp = -0.43 + 9.69 * (ctpoten**0.1875) 373c325,336 < endif ! if iamie=0 or not --- > if (ifcalccloc .eq. 0) then > theta0(isouth) = (-3.80+8.48*(ctpoten**0.1875))*dtr > theta0(inorth) = theta0(isouth) > offa(isouth) = 1.0*dtr > offa(inorth) = 1.0*dtr > dskofa(isouth) = 0. > dskofa(inorth) = 0. > arad(isouth) = max(rcp,rhp) > arad(inorth) = max(rcp,rhp) > rrad(isouth) = arad(isouth)*dtr > rrad(inorth) = arad(inorth)*dtr > endif 441,451c404,413 < if (iamie .eq. 0) then < offc(isouth) = 1.1*dtr < offc(inorth) = 1.1*dtr < dskofc(isouth) = (-0.08 + 0.15*byloc)*dtr < dskofc(inorth) = (-0.08 - 0.15*byloc)*dtr < phid(isouth) = (9.39 + 0.21*byloc - 12.) * h2deg * dtr < phid(inorth) = (9.39 - 0.21*byloc - 12.) * h2deg * dtr < phin(isouth) = (23.50 + 0.15*byloc - 12.) * h2deg * dtr < phin(inorth) = (23.50 - 0.15*byloc - 12.) * h2deg * dtr < endif ! iamie==0 < --- > if (ifcalccloc .eq. 0) then > offc(isouth) = 1.1*dtr > offc(inorth) = 1.1*dtr > dskofc(isouth) = (-0.08 + 0.15*byloc)*dtr > dskofc(inorth) = (-0.08 - 0.15*byloc)*dtr > phid(isouth) = (9.39 + 0.21*byloc - 12.) * h2deg * dtr > phid(inorth) = (9.39 - 0.21*byloc - 12.) * h2deg * dtr > phin(isouth) = (23.50 + 0.15*byloc - 12.) * h2deg * dtr > phin(inorth) = (23.50 - 0.15*byloc - 12.) * h2deg * dtr > endif 462,463d423 < < ! 467,468c427 < write(6,"(' iamie=',i2)") iamie < write(6,"(' cusp: alfac=',f8.3,' ec=',f8.3,' fc=',e10.4)") --- > write(6,"(' cusp: alfac=',f8.3,' ec=',f8.3,' fc=',e12.4)") 470c429 < write(6,"(' drizzle: alfad=',f8.3,' ed=',f8.3,' fd=',e10.4)") --- > write(6,"(' drizzle: alfad=',f8.3,' ed=',f8.3,' fd=',e12.4)") 477,478c436,437 < if (add_sproton) write(6,"(' solar protons: alfa_sp=',f8.3, < | ' e_sp = ',e10.4,' flx_sp = ',e10.4)") alfa_sp,e_sp,flx_sp --- > if (add_sproton) write(6,"(' solar protons: alfa_sp=',f10.3, > | ' e_sp = ',e10.3,' flx_sp = ',e10.3)") alfa_sp,e_sp,flx_sp 480,481c439,440 < | write(6,"(' high-energy electrons: alfa30 = ',f8.3, < | ' e30 = ',f8.3)") alfa30,e30 --- > | write(6,"(' high-energy electrons: alfa30 = ',f10.3, > | ' e30 = ',f10.3)") alfa30,e30 618,619d576 < ! (if this is an amie run, then theta0 and phid have been taken < ! from amie, see aurora_cons) 683,684c640 < real,parameter :: s10=0.174532925 < real :: clat --- > 748,770d703 < ! < ! Recalculate alfa, flux, drizl if an amie run (as in old sub amiepa): < if (iamie > 0) then < ! < ! Insure ekvg mean energy >= 1.: < ekvg(i,lat) = max(ekvg(i,lat),1.) < alfa(i+2) = ekvg(i,lat)/2. < flux(i+2) = efxg(i,lat)/(2.*alfa(i+2)*1.602e-9) < eflux(i+2) = flux(i+2)*fac_p2e*alfa(i+2) < clat = acos(sin(abs(dlat_aur(i)))) < alfa2(i+2) = alfa20 < ! drizl(i+2) = exp(-((clat-crad(ihem)+ < ! | abs(clat-crad(ihem)))/(2.*s10))**2) < ! flux2(i+2) = e20*(1.-re2*coslamda(i))* < ! | exp(-((clat-crad(ihem))/ < ! | halfwidth(i))**2) / (2.*alfa2(i+2)*1.602E-9) < drizl(i+2) = exp(-((clat-rradp(ihem)+ < | abs(clat-rradp(ihem)))/(1*s10))**2) < flux2(i+2) = e20*(1.-re2*coslamda(i))* < | exp(-((clat-rrad(ihem))/ < | halfwidth(i))**2) / (2.*alfa2(i+2)*1.602E-9) < eflux2(i+2) = flux2(i+2)*fac_p2e*alfa2(i+2) < endif 805a739 > ! NOTE: geng is divided by 2 from MIX before entering TIEGCM so this geng=characteristic 1055,1088c989,998 < ! Changed by G. Lu - 11/11/2011 < if (1.5*qo2p_aur(lev0,i) > 0.5*qo2p_aur(lev0+1,i)) then < qo2p(lev0,i,lat) = qo2p(lev0,i,lat)+1.5*qo2p_aur(lev0,i)- < | 0.5*qo2p_aur(lev0+1,i) < ! else < ! qo2p(lev0,i,lat)=qo2p(lev0,i,lat)+qo2p_aur(lev0,i) < endif < if (1.5*qop_aur (lev0,i) > 0.5*qop_aur (lev0+1,i)) then < qop (lev0,i,lat) = qop(lev0,i,lat)+1.5*qop_aur(lev0,i)- < | 0.5*qop_aur(lev0+1,i) < ! else < ! qop (lev0,i,lat) = qop(lev0,i,lat)+qop_aur(lev0,i) < endif < if (1.5*qn2p_aur(lev0,i)>0.5*qn2p_aur(lev0+1,i)) then < qn2p(lev0,i,lat) = qn2p(lev0,i,lat)+1.5*qn2p_aur(lev0,i)- < | 0.5*qn2p_aur(lev0+1,i) < ! else < ! qn2p(lev0,i,lat) = qn2p(lev0,i,lat)+qn2p_aur(lev0,i) < endif < if (1.5*qn2p_aur(lev0,i)>0.5*qn2p_aur(lev0+1,i)) then < qnp (lev0,i,lat) = qnp (lev0,i,lat)+ .22/.7 * < | (1.5*qn2p_aur(lev0,i)-0.5*qn2p_aur(lev0+1,i)) < ! else < ! qnp (lev0,i,lat) = qnp (lev0,i,lat)+.22/.7*qn2p_aur(lev0,i) < endif < if (1.5*qn2p_aur(lev0,i)>0.5*qn2p_aur(lev0+1,i)) then < qtef(lev0,i,lat)=qtef(lev0,i,lat)+1.57*(1.5*qn2p_aur(lev0,i)- < | 0.5*qn2p_aur(lev0+1,i)) < ! else < ! qtef(lev0,i,lat) = qtef(lev0,i,lat)+1.57*qn2p_aur(lev0,i) < endif < if (qn2p(lev0,i,lat) .lt. 0.) write(6,"('AURORA: ', < | 'i,lev0,lat,qn2p = ',3i3,4e12.4)") qo2p(lev0,i,lat), < | qop(lev0,i,lat),qn2p(lev0,i,lat),qnp (lev0,i,lat) --- > qo2p(lev0,i,lat) = qo2p(lev0,i,lat)+1.5*qo2p_aur(lev0,i)- > | 0.5*qo2p_aur(lev0+1,i) > qop (lev0,i,lat) = qop (lev0,i,lat)+1.5*qop_aur (lev0,i)- > | 0.5*qop_aur (lev0+1,i) > qn2p(lev0,i,lat) = qn2p(lev0,i,lat)+1.5*qn2p_aur(lev0,i)- > | 0.5*qn2p_aur(lev0+1,i) > qnp (lev0,i,lat) = qnp (lev0,i,lat)+ .22/.7 * > | (1.5*qn2p_aur(lev0,i)-0.5*qn2p_aur(lev0+1,i)) > qtef(lev0,i,lat) = qtef(lev0,i,lat)+1.57*(1.5*qn2p_aur(lev0,i)- > | 0.5*qn2p_aur(lev0+1,i))