netcdf dB_t_M01 {
dimensions:
	time = UNLIMITED ; // (27 currently)
	mlon1 = 100 ;
	mlon2 = 100 ;
	mlat1 = 161 ;
	mlat2 = 160 ;
	mlat3qd = 162 ;
	hgt = 54 ;
	hgt_r = 55 ;
	lon_qd = 100 ;
	lat_qd = 220 ;
	glon = 73 ;
	glat = 91 ;
	ghgt = 55 ;
	fldpts_qdlat1 = 5886 ;
	fldpts_hgt1 = 5886 ;
	fldpts_qdlat2 = 5778 ;
	fldpts_hgt2 = 5778 ;
	fldpts_qdlat3 = 5940 ;
	fldpts_hgt3 = 5940 ;
	fldpts_qdlatp = 5886 ;
	fldpts_hgtp = 5886 ;
variables:
	double time(time) ;
	double mlon1(mlon1) ;
	double mlon2(mlon2) ;
	double lon_qd(lon_qd) ;
	double mlat1(mlat1) ;
	double mlat2(mlat2) ;
	double lat_qd(lat_qd) ;
	double mlat3qd(mlat3qd) ;
	double hgt(hgt) ;
	double hgt_r(hgt_r) ;
	double fldpts_qdlat1(fldpts_qdlat1) ;
	double fldpts_hgt1(fldpts_hgt1) ;
	double fldpts_qdlat2(fldpts_qdlat2) ;
	double fldpts_hgt2(fldpts_hgt2) ;
	double fldpts_qdlat3(fldpts_qdlat3) ;
	double fldpts_hgt3(fldpts_hgt3) ;
	double fldpts_qdlatp(fldpts_qdlatp) ;
	double fldpts_hgtp(fldpts_hgtp) ;
	double glon(glon) ;
	double glat(glat) ;
	double ghgt(ghgt) ;
	
	
	double year(time) ;
		year:long_name = "year" ;
	double doy(time) ;
		doy:long_name = "doy" ;
	double f107(time) ;
		f107:long_name = "f107" ;
	double ap(time) ;
		ap:long_name = "ap" ;
	double ut(time) ;
		ut:long_name = "ut" ;
		ut:units = "hours" ;
	double ctpoten(time) ;
		ctpoten:long_name = "ctpo" ;
	double hpower(time) ;
		hpower:long_name = "hpow" ;

! Jf1(i-0.5,l,k) = 0.5*I1qd(i-0.5,l,k)/r(k)^1.5/(r(k+0.5)^0.5-r(k-0.5)^0.5)/dlamq
! Jf2(i,l+0.5,k) = I2qd(i,l+0.5,k)/r(k)/(r(k+0.5)-r(k-0.5))/dlonq/cos(latq(l+0.5))
! Jr(i,l,k-0.5)  = I3qd(i,l,k-0.5)/M3qd(i,l,k-0.5)		
	double Jf1(time, hgt, lat_qd, lon_qd) ;
		Jf1:long_name = "Jf1" ;
		Jf1:units = "A/m2" ;
	double Jf2(time, hgt, lat_qd, lon_qd) ;
		Jf2:long_name = "Jf2" ;
		Jf2:units = "A/m2" ;
	double Jr(time, hgt, lat_qd, lon_qd) ;
		Jr:long_name = "Jr " ;
		Jr:units = "A/m2" ;
		
! calculate Jeej = (Jf1hor*f1+Jf1hor*f2)*f1/|f1|
	double Jeej(time, hgt, lat_qd, lon_qd) ;
		Jeej:long_name = "Jeej  " ;
		Jeej:units = "A/m2" ;
	    
! calculate Jfac = Jf1_hor*f1*bo^+Jf1_hor*f2*bo^+ Jr*k*bo^
! F1 and F2 are horizontal since g3=Fk and f1=g2 x g3 and f2=g3 x g1		
	double Jfac(time, hgt, lat_qd, lon_qd) ;
		Jfac:long_name = "Jfac  " ;
		Jfac:units = "A/m2" ;

! calculate current Je2 Je2 = J*d2 = (Jf1hor*f1+Jf2hor*f2+Jr*k)*d2		
	double Je2J(time, hgt, lat_qd, lon_qd) ;
		Je2J:long_name = "Je2J  " ;
		Je2J:units = "A/m2" ;

! alterative calculation: Je2 = -(R/r)cosI_m Jr - (r/R)^1.5*D*sinI*Jf2		 
	double Je2JA(time, hgt, lat_qd, lon_qd) ;
		Je2JA:long_name = "Je2JA  " ;
		Je2JA:units = "A/m2" ;
! alterative calculation: Je2 = -(R/r)cosI_m Jr 
	double Je2JA1(time, hgt, lat_qd, lon_qd) ;
		Je2JA1:long_name = "Je2JA1" ;
		Je2JA1:units = "A/m2" ;
! alterative calculation: Je2 = - (r/R)^1.5*D*sinI*Jf2		
	double Je2JA2(time, hgt, lat_qd, lon_qd) ;
		Je2JA2:long_name = "Je2JA2" ;
		Je2JA2:units = "A/m2" ;

! calculate horizontal current J = Jf1hor*f1+Jf2hor*f2+Jr*k 
! Jf1hor = Jf1-(k^ dot g1)Jr (eq. 261)
! Jf2hor = Jf2-(k^ dot g2)Jr (eq. 262)		
	double Jf1hor(time, hgt, lat_qd, lon_qd) ;
		Jf1hor:long_name = "Jf1hor" ;
		Jf1hor:units = "A/m2" ;
	double Jf2hor(time, hgt, lat_qd, lon_qd) ;
		Jf2hor:long_name = "Jf2" ;
		Jf2hor:units = "A/m2" ;

! mapped fiels		
	double un_s1(time, fldpts_qdlat1, mlon1) ;
		un_s1:long_name = "un_s1" ;
		un_s1:units = "m/s" ;
	double un_s2(time, fldpts_qdlat2, mlon2) ;
		un_s2:long_name = "un_s2" ;
		un_s2:units = "m/s" ;
	double vn_s1(time, fldpts_qdlat1, mlon1) ;
		vn_s1:long_name = "vn_s1" ;
		vn_s1:units = "m/s" ;
	double vn_s2(time, fldpts_qdlat2, mlon2) ;
		vn_s2:long_name = "vn_s2" ;
		vn_s2:units = "m/s" ;
	double sigH_s1(time, fldpts_qdlat1, mlon1) ;
		sigH_s1:long_name = "sigH_s1" ;
		sigH_s1:units = "S/m" ;
	double sigP_s1(time, fldpts_qdlat1, mlon1) ;
		sigP_s1:long_name = "sigP_s1" ;
		sigP_s1:units = "S/m" ;
	double sigP_s2(time, fldpts_qdlat2, mlon2) ;
		sigP_s2:long_name = "sigP_s2" ;
		sigP_s2:units = "S/m" ;
	double sigH_s2(time, fldpts_qdlat2, mlon2) ;
		sigH_s2:long_name = "sigH_s2" ;
		sigH_s2:units = "S/m" ;
	double J3lb(time, mlat3qd, mlon2) ;
		J3lb:long_name = "J3lb" ;
		J3lb:units = "A/m2" ;
		
! sinI inclination		
	double sinI_s1(time, fldpts_qdlat1, mlon1) ;
		sinI_s1:long_name = "sinI_s1" ;
		sinI_s1:units = "[-]" ;
	double sinI_s2(time, fldpts_qdlat2, mlon2) ;
		sinI_s2:long_name = "sinI_s2" ;
		sinI_s2:units = "[-]" ;
		
! M and N values needed for solving and later calculation of J		
	double N1p(time, fldpts_qdlat1, mlon1) ;
		N1p:long_name = "N1p" ;
		N1p:units = "?" ;
	double N1h(time, fldpts_qdlat1, mlon1) ;
		N1h:long_name = "N1h" ;
		N1h:units = "?" ;
	double M1(time, fldpts_qdlat1, mlon1) ;
		M1:long_name = "M1" ;
		M1:units = "?" ;
	double N2p(time, fldpts_qdlat2, mlon2) ;
		N2p:long_name = "N2p" ;
		N2p:units = "?" ;
	double N2h(time, fldpts_qdlat2, mlon2) ;
		N2h:long_name = "N2h" ;
		N2h:units = "?" ;
	double M2(time, fldpts_qdlat2, mlon2) ;
		M2:long_name = "M2" ;
		M2:units = "?" ;
	double M3(time, fldpts_qdlat3, mlon2) ;
		M3:long_name = "M3" ;
		M3:units = "?" ;
		
	double Ne1(time, fldpts_qdlat1, mlon1) ;
		Ne1:long_name = "Ne1" ;
		Ne1:units = "#/m3" ;
	double Tei1(time, fldpts_qdlat1, mlon1) ;
		Tei1:long_name = "Tei1" ;
		Tei1:units = "K   " ;
	double Ne2(time, fldpts_qdlat2, mlon2) ;
		Ne2:long_name = "Ne2" ;
		Ne2:units = "#/m3" ;
	double Tei2(time, fldpts_qdlat2, mlon2) ;
		Tei2:long_name = "Tei2" ;
		Tei2:units = "K   " ;

!
! equation (83) page 10 Art's script
! I1(i+0.5,j,k) = N1p(i+0.5,j,k)*[Phi(i,j)-Phi(i+1,j)]
!      -N1h(i+0.5,j,k)[Phi(i,j-1)+Phi(i+1,j-1)-Phi(i,j+1)-Phi(i+1,j+1)]
!      +M1(i+0.5,j,k)Je1D(i+0.5,j,k)
!
! equation (85) page 10 Art's script
! I2(i,j+0.5,k) = N2h(i,j+0.5,k)*[Phi(i-1,j)+Phi(i-1,j+1)-Phi(i+1,j)-Phi(i+1,j+1)]
!      +N2p(i,j+0.5,k)[Phi(i,j)-Phi(i,j+1)]
!      +M2(i,j+0.5,k)Je2D(i,j+0.5,k)
!
! equation (63') page 7 Art's script
! for all i, and j=2,nmlat_h
!      I3(i,j,k+0.5) = I3(i,j,k-0.5) + I1(i-0.5,j,k)-I1(i+0.5,j,k)+I2(i,j-0.5,k)-I2(i,j+0.5,k)		
	double I1(time, fldpts_qdlat1, mlon1) ;
		I1:long_name = "I1" ;
		I1:units = "A" ;
	double I1_1(time, fldpts_qdlat1, mlon1) ;
		I1_1:long_name = "I1_1" ;
		I1_1:units = "A" ;
	double I1_2(time, fldpts_qdlat1, mlon1) ;
		I1_2:long_name = "I1_2" ;
		I1_2:units = "A" ;
	double I1_3(time, fldpts_qdlat1, mlon1) ;
		I1_3:long_name = "I1_3" ;
		I1_3:units = "A" ;
	double I2(time, fldpts_qdlat2, mlon2) ;
		I2:long_name = "I2" ;
		I2:units = "A" ;
	double I2_1(time, fldpts_qdlat2, mlon2) ;
		I2_1:long_name = "I2_1" ;
		I2_1:units = "A" ;
	double I2_2(time, fldpts_qdlat2, mlon2) ;
		I2_2:long_name = "I2_2" ;
		I2_2:units = "A" ;
	double I2_3(time, fldpts_qdlat2, mlon2) ;
		I2_3:long_name = "I2_3" ;
		I2_3:units = "A" ;
	double I2oM2(time, fldpts_qdlat2, mlon2) ;
		I2oM2:long_name = "I2oM2" ;
		I2oM2:units = "A/m2" ;
	double I3(time, fldpts_qdlat3, mlon2) ;
		I3:long_name = "I3" ;
		I3:units = "A" ;
		
S is the wind driven and ionospheric current sources (89) 
! S(i,j,k) = M1(i-0.5,j,k)*Je1D(i-0.5,j,k)-M1(i+0.5,j,k)*Je1D(i+0.5,j,k) +
!       M2(i,j-0.5,k)*J2D(i,j-0.5,k)-M2(i,j+0.5,k)*Je1D(i,j+0.5,k)		
	double S(time, fldpts_qdlatp, mlon1) ;
		S:long_name = "S" ;
		S:units = "?" ;

! diagnostic to get wind driven current in f coordinate system
! process for S2 points
! 1. calculate J|| to make local current balance no vertical current (with Je3 e3= J|| b^)
!   ( Je1D e1 + Je2D e2+ J||b^) dot k^ = 0 -> J||= -(Je1D e1 dot k^ + Je2D e2 dot k^)/b^ dot k^
! 2. calcluate Jf1D and Jf2D and JrD  
!   Jf1 = g1 dot(Je1D e1 + Je2D e2+ J||b^)  horizontal component
!   Jf2 = g2 dot(Je1D e1 + Je2D e2+ J||b^)  horizontal component       		
	double Jf1Dyn(time, fldpts_qdlat2, mlon2) ;
		Jf1Dyn:long_name = "Jf1Dyn" ;
		Jf1Dyn:units = "A/m2" ;
	double Jf2Dyn(time, fldpts_qdlat2, mlon2) ;
		Jf2Dyn:long_name = "Jf2Dyn" ;
		Jf2Dyn:units = "A/m2" ;

! diagnostic to get wind driven current 		
	double Je1D(time, fldpts_qdlat1, mlon1) ;
		Je1D:long_name = "Je1D" ;
		Je1D:units = "?" ;
	double Je2D(time, fldpts_qdlat2, mlon2) ;
		Je2D:long_name = "Je2D" ;
		Je2D:units = "?" ;
		
!   Jf1(Je2p) = g1 dot(Je2p e2)  horizontal component
!   Jf2(Je2p) = g2 dot(Je2p e2)  horizontal component
	double Jf1Ion2(time, fldpts_qdlat2, mlon2) ;
		Jf1Ion2:long_name = "Jf1Ion2" ;
		Jf1Ion2:units = "A/m2" ;
	double Jf2Ion2(time, fldpts_qdlat2, mlon2) ;
		Jf2Ion2:long_name = "Jf2Ion2" ;
		Jf2Ion2:units = "A/m2" ;
		
	double Je1Ion(time, fldpts_qdlat1, mlon1) ;
		Je1Ion:long_name = "Je1Ion" ;
		Je1Ion:units = "A/m2" ;
	double Je2Ion(time, fldpts_qdlat2, mlon2) ;
		Je2Ion:long_name = "Je2Ion" ;
		Je2Ion:units = "A/m2" ;
		
! electric potential		
	double poten(time, mlat1, mlon2) ;
		poten:long_name = "poten" ;
		poten:units = "V" ;
! field-aligned current density based on high latitude prescribed potential		
	double fac_dens(time, mlat1, mlon2) ;
		fac_dens:long_name = "fac_dens" ;
		fac_dens:units = "A/m2" ;
! high latitude electric potential
	double poten_hl(time, mlat1, mlon1) ;
		poten_hl:long_name = "poten_hl" ;
		poten_hl:units = "V" ;

! conductances		
	double ZigP1(time, mlat1, mlon1) ;
		ZigP1:long_name = "ZigP1" ;
		ZigP1:units = "S" ;
	double ZigP2(time, mlat2, mlon2) ;
		ZigP2:long_name = "ZigP2" ;
		ZigP2:units = "S" ;
	double ZigH1(time, mlat1, mlon1) ;
		ZigH1:long_name = "ZigH1" ;
		ZigH1:units = "S" ;
	double ZigH2(time, mlat2, mlon2) ;
		ZigH2:long_name = "ZigH2" ;
		ZigH2:units = "S" ;
		
! electric field		
	double Ed1_S1(time, mlat1, mlon1) ;
		Ed1_S1:long_name = "Ed1_S1" ;
		Ed1_S1:units = "V/m" ;
	double Ed2_S1(time, mlat1, mlon1) ;
		Ed2_S1:long_name = "Ed2_S1" ;
		Ed2_S1:units = "V/m" ;
	double Ed1_S2(time, mlat2, mlon2) ;
		Ed1_S2:long_name = "Ed1_S2" ;
		Ed1_S2:units = "V/m" ;
	double Ed2_S2(time, mlat2, mlon2) ;
		Ed2_S2:long_name = "Ed2_S2" ;
		Ed2_S2:units = "V/m" ;

! ExB drift		
	double Ve1_S1(time, mlat1, mlon1) ;
		Ve1_S1:long_name = "Ve1_S1" ;
		Ve1_S1:units = "m/s" ;
	double Ve2_S1(time, mlat1, mlon1) ;
		Ve2_S1:long_name = "Ve2_S1" ;
		Ve2_S1:units = "m/s" ;

! dB calculation		
	double Be_grd(time, glat, glon) ;
		Be_grd:long_name = "Be_grd" ;
		Be_grd:units = "T" ;
	double Bn_grd(time, glat, glon) ;
		Bn_grd:long_name = "Bn_grd" ;
		Bn_grd:units = "T" ;
	double Bu_grd(time, glat, glon) ;
		Bu_grd:long_name = "Bu_grd" ;
		Bu_grd:units = "T" ;
	double Be_LEO(time, glat, glon) ;
		Be_LEO:long_name = "Be_LEO" ;
		Be_LEO:units = "T" ;
	double Bn_LEO(time, glat, glon) ;
		Bn_LEO:long_name = "Bn_LEO" ;
		Bn_LEO:units = "T" ;
	double Bu_LEO(time, glat, glon) ;
		Bu_LEO:long_name = "Bu_LEO" ;
		Bu_LEO:units = "T" ;
	double Be_70W(time, ghgt, glat) ;
		Be_70W:long_name = "Be_70W" ;
		Be_70W:units = "T" ;
	double Bn_70W(time, ghgt, glat) ;
		Bn_70W:long_name = "Bn_70W" ;
		Bn_70W:units = "T" ;
	double Bu_70W(time, ghgt, glat) ;
		Bu_70W:long_name = "Bu_70W" ;
		Bu_70W:units = "T" ;

! equivalent current function		
	double Psi_geo(time, glat, glon) ;
		Psi_geo:long_name = "Psi_geo" ;
		Psi_geo:units = "A" ;
}
