#include "dims.h" ! gl 23/07/2002 save extra fields into secondary history ! subroutine extraoutput use cons_module,only: pi, dlamda, dphi, sn implicit none ! #include "params.h" #include "index.h" #include "buff.h" #include "phys.h" #include "dynamo.h" #include "trgm.h" #include "extrafields.h" ! ! Local: integer :: j1,i,lonm,lonp,k,lat real :: re CCCCCCCCCCCCCCCC calculate the field-aligned currents CCCCCCCCCCCCCCCCCC re = 6.481e+6 fwind(:,:) = 0. famie(:,:) = 0. ftot(:,:) = 0. do 5150 lat = 2, zjmx-1 j1 = zjmx+1-lat if (abs(sin(dipmag(i,lat))) .lt. 0.2) then c write(6,"('extroutput: small magdip angle at ', c + 'glat,sndip(i,lat)=',f6.2,f7.4)") c + -92.5+dphi*lat,sin(dipmag(i,lat)) go to 5150 endif do i=1, zimx+1 lonm = i -1 if (lonm .eq. 0) lonm = zimx+1 lonp = i + 1 if (lonp .eq. zimx+2) lonp = 1 C Field-aligned current from wind fwind(i,lat) = ((fwindu(lonp,lat)-fwindu(lonm,lat))/ 1 dlamda + (sn(lat+1)*fwindv(i,lat+1)-sn(lat-1)* 2 fwindv(i,lat-1))/dphi) / (2.*re*sn(lat)) / 3 (-sin(dipmag(i,lat))) C Field-aligned current from ion drift famie(i,lat) = ((famieu(lonp,lat)-famieu(lonm,lat))/ 1 dlamda + (sn(lat+1)*famiev(i,lat+1)-sn(lat-1)* 2 famiev(i,lat-1))/dphi) / (2.*re*sn(lat)) / 3 (-sin(dipmag(i,lat))) C Total field-aligned currents ftot(i,lat) = famie(i,lat) + fwind(i,lat) C *** calculate height-profiles of field-aligned currents *** do k=1,zkmx fwind_sec(i,lat,k) = + ((fwindu_sec(lonp,lat,k)-fwindu_sec(lonm,lat,k))/ + dlamda + (sn(lat+1)*fwindv_sec(i,lat+1,k)-sn(lat-1)* + fwindv_sec(i,lat-1,k))/dphi) / (2.*re*sn(lat)) / + (-sin(dipmag(i,lat))) famie_sec(i,lat,k) = + ((famieu_sec(lonp,lat,k)-famieu_sec(lonm,lat,k))/ + dlamda + (sn(lat+1)*famiev_sec(i,lat+1,k)-sn(lat-1)* + famiev_sec(i,lat-1,k))/dphi) / (2.*re*sn(lat)) / + (-sin(dipmag(i,lat))) ftot_sec(i,lat,k) = famie_sec(i,lat,k)+fwind_sec(i,lat,k) enddo C *** end of calculate height-profiles of field-aligned currents *** enddo famie_sec(:,lat,zkmxp) = famie(:,lat) fwind_sec(:,lat,zkmxp) = fwind(:,lat) 5150 continue fwind(:,1) = (fwind(:,2)-fwind(:,1)) / re/dphi fwind(:,zjmx) = (fwind(:,zjmx)-fwind(:,zjmx-1))/re/dphi famie(:,1) = (famie(:,2)-famie(:,1)) / re/dphi famie(:,zjmx) = (famie(:,zjmx)-famie(:,zjmx-1))/re/dphi ftot(:,1) = fwind(:,1) + famie(:,1) ftot(:,zjmx) = fwind(:,zjmx) + famie(:,zjmx) do k=1,zkmx fwind_sec(:,1,k)=(fwind_sec(:,2,k)-fwind_sec(:,1,k))/re/dphi fwind_sec(:,zjmx,k) = (fwind_sec(:,zjmx,k)- + fwind_sec(:,zjmx-1,k))/re/dphi famie_sec(:,1,k)=(famie_sec(:,2,k)-famie_sec(:,1,k))/re/dphi famie_sec(:,zjmx,k) = (famie_sec(:,zjmx,k)- + famie_sec(:,zjmx-1,k))/re/dphi ftot_sec(:,1,k) = fwind_sec(:,1,k)+famie_sec(:,1,k) ftot_sec(:,zjmx,k)=fwind_sec(:,zjmx,k)+famie_sec(:,zjmx,k) enddo famie_sec(:,1,zkmxp) = famie(:,1) fwind_sec(:,1,zkmxp) = fwind(:,1) famie_sec(:,zjmx,zkmxp) = famie(:,zjmx) fwind_sec(:,zjmx,zkmxp) = fwind(:,zjmx) CCCCCCCCCCCCCCCCCCC end of FAC calculation CCCCCCCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCC save secondary histories CCCCCCCCCCCCCCCCCCCCC ! ! Add sigmas to secondary histories: ! subroutine addfsech(fname,f,idim1,idim2,ndim2,ilat) ! Field passed to addfsech must have ZIMXP (76) in 1st dimension, ! and ZKMXP in 2nd dimension. ! Since sigma1,2 have zimxp1 (73) in 1st dimension, so transfer ! to local arrays sigma1_sech and sigma2_sech w/ periodic points. ! Since ndim2 < idim2 by 1 (ZKMX < ZKMXP) in this case, addfsech will ! pad the top level with spval. ! ZIMXP = 76, ZIMXP1 = 72+1 = 73, ZKMX = 28, ZKMXP = ZKMX+1 = 29 ! From dynamo.h: SIGMA1(ZIMXP1,ZJMX,ZKMX),SIGMA2(ZIMXP1,ZJMX,ZKMX) ! do lat = 1, zjmx sigma1_sech(1:zimxp1,1:zkmx) = sigma1(:,lat,:) sigma1_sech(1:zimxp1,zkmxp) = sigp(:,lat) sigma1_sech(zimxp1+1,:) = sigma1_sech(2,:) ! 74 <- 2 sigma1_sech(zimxp1+2,:) = sigma1_sech(3,:) ! 75 <- 3 sigma1_sech(zimxp1+3,:) = sigma1_sech(4,:) ! 76 <- 4 call addfsech('PEDERSEN',' ',' ',sigma1_sech,zimxp, | zkmxp,zkmxp,lat) sigma2_sech(1:zimxp1,1:zkmx) = sigma2(:,lat,:) sigma2_sech(1:zimxp1,zkmxp) = sigh(:,lat) sigma2_sech(zimxp1+1,:) = sigma2_sech(2,:) ! 74 <- 2 sigma2_sech(zimxp1+2,:) = sigma2_sech(3,:) ! 75 <- 3 sigma2_sech(zimxp1+3,:) = sigma2_sech(4,:) ! 76 <- 4 call addfsech('HALL',' ',' ',sigma2_sech,zimxp,zkmxp,zkmxp,lat) qwind_sech(1:zimxp1,1:zkmxp) = qwind_sec(:,lat,:) qwind_sech(zimxp1+1,:) = qwind_sech(2,:) ! 74 <- 2 qwind_sech(zimxp1+2,:) = qwind_sech(3,:) ! 75 <- 3 qwind_sech(zimxp1+3,:) = qwind_sech(4,:) ! 76 <- 4 call addfsech('QWIND',' ',' ',qwind_sech,zimxp,zkmxp,zkmxp,lat) qamie_sech(1:zimxp1,1:zkmxp) = qamie_sec(:,lat,:) qamie_sech(zimxp1+1,:) = qamie_sech(2,:) ! 74 <- 2 qamie_sech(zimxp1+2,:) = qamie_sech(3,:) ! 75 <- 3 qamie_sech(zimxp1+3,:) = qamie_sech(4,:) ! 76 <- 4 call addfsech('QAMIE',' ',' ',qamie_sech,zimxp,zkmxp,zkmxp,lat) fwind_sech(1:zimxp1,1:zkmxp) = fwind_sec(:,lat,:) fwind_sech(zimxp1+1,:) = fwind_sech(2,:) ! 74 <- 2 fwind_sech(zimxp1+2,:) = fwind_sech(3,:) ! 75 <- 3 fwind_sech(zimxp1+3,:) = fwind_sech(4,:) ! 76 <- 4 call addfsech('FWIND',' ',' ',fwind_sech,zimxp,zkmxp,zkmxp,lat) famie_sech(1:zimxp1,1:zkmxp) = famie_sec(:,lat,:) famie_sech(zimxp1+1,:) = famie_sech(2,:) ! 74 <- 2 famie_sech(zimxp1+2,:) = famie_sech(3,:) ! 75 <- 3 famie_sech(zimxp1+3,:) = famie_sech(4,:) ! 76 <- 4 call addfsech('FAMIE',' ',' ',famie_sech,zimxp,zkmxp,zkmxp,lat) work_sech(1:zimxp1,1:zkmxp) = work_sec(:,lat,:) work_sech(zimxp1+1,:) = work_sech(2,:) ! 74 <- 2 work_sech(zimxp1+2,:) = work_sech(3,:) ! 75 <- 3 work_sech(zimxp1+3,:) = work_sech(4,:) ! 76 <- 4 call addfsech('WORK',' ',' ',work_sech,zimxp,zkmxp,zkmxp,lat) wtot_sech(1:zimxp1,1:zkmxp) = wtot_sec(:,lat,:) wtot_sech(zimxp1+1,:) = wtot_sech(2,:) ! 74 <- 2 wtot_sech(zimxp1+2,:) = wtot_sech(3,:) ! 75 <- 3 wtot_sech(zimxp1+3,:) = wtot_sech(4,:) ! 76 <- 4 call addfsech('WTOT',' ',' ',wtot_sech,zimxp,zkmxp,zkmxp,lat) tec_sech(1:zimxp1,1:zkmxp) = tec_sec(:,lat,:) tec_sech(zimxp1+1,:) = tec_sech(2,:) ! 74 <- 2 tec_sech(zimxp1+2,:) = tec_sech(3,:) ! 75 <- 3 tec_sech(zimxp1+3,:) = tec_sech(4,:) ! 76 <- 4 call addfsech('TEC',' ',' ',tec_sech,zimxp,zkmxp,zkmxp,lat) enddo CCCCCCCCCCCCCCCCCCC end of save secondary histories CCCCCCCCCCCCCCCCCCCCC end subroutine extraoutput