SUBROUTINE tomhd (rec, ierr) USE Rcm_mod_subs, ONLY : isize, jsize, ksize, kcsize,jwrap, nptmax, & colat, aloct, itrack, v, birk, & bmin, xmin, ymin, vm, pmin, rmin, & eeta, alamc, etac, ikflavc, fudgec, & bi, bj, etab, ikflav, mpoint, npoint, & birk_avg, v_avg, eeta_avg, & pi, read_array, label, LUN, & boundary, bndloc, pressrcm,& Read_grid, Read_plasma,Get_boundary, & xmass, densrcm,imin_j,iprec,rprec USE intermediate_grid USE lfm_transfergrid USE intercomm_module, ONLY : ExportMagnetosphere USE constants, ONLY : mass_proton USE rice_housekeeping_module ! !============================================================== ! purpose: ! To read rcm data at a specified record (rec) and place ! into a common block '/press2_eq/' 2d pressure information ! now replaced - frt 10.03 ! It also writes to files: ! 'rcmu.dat' = unformatted time-averaged rcm data for plotting ! 'press.dat' = formatted pressure data in the equat plane for plots ! ! inputs: ! rec = record to read rcm data ! ! ! program to read rcm output files and write out results in ! this version also spits out pressure contours in the equatorial ! plane for input to setupdg.f which setups the pressure for the ! multi channel version 02/00 frt ! ! 4/18/95 rws ! 9/18/95 frt ! 9/1/98 ythis version takes into account the the rcm record ! can differ to the record ! 7/09 -restructured to use modules and allow to transfer ! the LFM grid - frt ! ! !============================================================== ! IMPLICIT NONE INTEGER(iprec), INTENT (IN) :: rec INTEGER(iprec), INTENT (OUT) :: ierr ! INTEGER(iprec) :: iunit, iunit1, iunit2, iout, ig, jg INTEGER(iprec) :: kkk, i, j, L, k, itime, n INTEGER(iprec), PARAMETER :: idim=isize, jdim=jsize, jdimr=jsize REAL(rprec), PARAMETER :: vmmax = 2000.0 REAL(rprec) :: bbi, bbj, di, dj, & emp, vmp, ptemp, rhotemp,ctemp REAL(rprec) :: xval,yval,pval,rval LOGICAL :: L_flag ! real(rprec), parameter:: pressmin0 =1.0e-12 ! real(rprec),parameter :: densitymin0 =0.1*1e6*1.67e-27 ! 0.1 ple/cc real(rprec) :: xg,yg real(rprec) :: xn,yn,zn real(rprec) :: xtrans,ytrans,ztrans,rtrans real(rprec) :: rinter real(rprec) :: peq,denseq real(rprec) :: max_xmin,min_xmin,max_ymin,min_ymin real(rprec), PARAMETER :: rtrans_min = 2.0 ! min radius to nudge back real(rprec), PARAMETER :: empmin = 1.0E+10 integer(iprec) :: in_tot,ierr2,mask integer(iprec) :: istart,iend,iit,jstart,jend,jit,kstart,kend,kit integer(iprec) :: ier,ifound integer(iprec) :: nxp,nyp integer(iprec) :: bnd_max integer(iprec), parameter :: i_offset = 1, imin_grid = 2 real(rprec) :: dens_plasmasphere logical, parameter :: use_plasmasphere = .true. logical, parameter :: use_ionosphere_interpolation = .true. save ig,jg,bbi,bbj,di,dj,iout,ierr2 save pval,rval,ptemp,rhotemp save xn,yn,zn,xtrans,ytrans,ztrans save ifound,mask INTEGER (iprec) :: itime0 INTEGER (iprec), ALLOCATABLE :: mpoint0(:), npoint0(:) REAL (rprec), ALLOCATABLE :: v0(:,:), birk0(:,:), vm0(:,:), bmin0(:,:),& xmin0(:,:), ymin0(:,:), rmin0(:,:), pmin0(:,:), eeta0(:,:,:),& bi0(:),bj0(:),etab0(:),v_avg0(:,:),birk_avg0(:,:),eeta_avg0(:,:,:) LOGICAL,PARAMETER :: avoid_boundaries = .false. INTEGER :: im,ip,jm,jp,km,kp INCLUDE 'rcmdir.h' ! Logic of this thing is such that rec must be >= 2 (RCM had to run at least once): IF (rec < 2) STOP ' tomhd called before RCM was called, aborting...' ! At this point, we assume that RCM arrays are all populated ! (plasma, b-field, bndloc, etc.): IF (L_write_rcmu) then IF (.NOT.allocated(v0) ) then ALLOCATE (v0(isize,jsize), birk0(isize,jsize), vm0(isize,jsize),& bmin0(isize,jsize),xmin0(isize,jsize), ymin0(isize,jsize),& rmin0(isize,jsize),pmin0(isize,jsize), & birk_avg0(isize,jsize), v_avg0(isize,jsize), & eeta0(isize,jsize,ksize), eeta_avg0(isize,jsize,ksize), & bi0(nptmax),bj0(nptmax),etab0(nptmax), & mpoint0(SIZE(mpoint)), npoint0(SIZE(npoint))) END IF ! save current RCM arrays so that another set can be retrieved then will put them back. v0 = v birk0 = birk vm0 = vm bmin0 = bmin xmin0 = xmin ymin0 = ymin rmin0 = rmin pmin0 = pmin eeta0 = eeta bi0 = bi bj0 = bj etab0 = etab mpoint0 = mpoint npoint0 = npoint v_avg0 = v_avg birk_avg0 = birk_avg eeta_avg0 = eeta_avg itime0 = label%intg(6) L = rec - 1 CALL Read_array (rcmdir//'rcmv' , L, label, ARRAY_2D = v) CALL Read_array (rcmdir//'rcmbirk', L, label, ARRAY_2D = birk) CALL Read_array (rcmdir//'rcmvm', L, label, ARRAY_2D = vm) CALL Read_array (rcmdir//'rcmbmin', L, label, ARRAY_2D = bmin) CALL Read_array (rcmdir//'rcmxmin', L, label, ARRAY_2D = xmin) CALL Read_array (rcmdir//'rcmymin', L, label, ARRAY_2D = ymin) rmin = SQRT (xmin**2+ymin**2) WHERE (xmin == 0.0 .AND. ymin == 0.0) pmin = 0.0 ELSEWHERE pmin = ATAN2 (ymin, xmin) END WHERE WHERE (pmin < 0.0) pmin = pmin + 2.0*pi CALL Read_array (rcmdir//'rcmeeta', L, label, ARRAY_3D = eeta) CALL Read_array (rcmdir//'rcmbi', L, label, ARRAY_1D = bi) CALL Read_array (rcmdir//'rcmbj', L, label, ARRAY_1D = bj) CALL Read_array (rcmdir//'rcmetab', L, label, ARRAY_1D = etab) CALL Read_array (rcmdir//'rcmitrack', L, label, ARRAY_1D = itrack) CALL Read_array (rcmdir//'rcmmpoint', L, label, ARRAY_1D = mpoint) CALL Read_array (rcmdir//'rcmnpoint', L, label, ARRAY_1D = npoint) CALL Read_array (rcmdir//'rcmvavg' , L, label, ARRAY_2D = v_avg) CALL Read_array (rcmdir//'rcmbirkavg', L, label, ARRAY_2D = birk_avg) CALL Read_array (rcmdir//'rcmeetaavg', L, label, ARRAY_3D = eeta_avg) itime = label%intg(6) WRITE (*,'(A,I9.9,A,I5.5)') 'TOMHD: Read RCM, T=',itime,', REC=',L IF (L_write_rcmu) call write_rcmu (L,0) v = v0 birk = birk0 vm = vm0 bmin = bmin0 xmin = xmin0 ymin = ymin0 rmin = rmin0 pmin = pmin0 eeta = eeta0 bi = bi0 bj = bj0 etab = etab0 mpoint = mpoint0 npoint = npoint0 v_avg = v_avg0 birk_avg = birk_avg0 eeta_avg = eeta_avg0 label%intg(6) = itime0 L = rec call write_rcmu (L,0) END IF ! Compute rcm pressure and density (on the ionospheric RCM grid): ! DO j = 1, jdim DO i = 1, idim pressrcm (i,j) = 0.0 densrcm (i,j) = 0.0 IF (vm(i,j) < 0.0) CYCLE DO k = 1, kcsize pressrcm(i,j) = pressrcm(i,j) + & 1.674E-35*ABS(alamc(k))*eeta(i,j,k)*vm(i,j)**2.5 ! in pascals ! normalize everything to the mass_proton, otherwise answer is below ! floating point minimum answer and gets zero in ples/m^3 densrcm(i,j) = densrcm(i,j) + & 1.5694E-16/mass_proton*xmass(ikflavc(k))*eeta(i,j,k)*vm(i,j)**1.5 END DO if(use_plasmasphere)then ! add a simple plasmasphere model based on carpenter 1992 or gallagher 2002 in ples/cc ! call carpenter(rmin(i,j),dens_plasmasphere) call gallagher(rmin(i,j),dens_plasmasphere) densrcm(i,j) = densrcm(i,j) + dens_plasmasphere*1.0e6 end if END DO END DO max_xmin = maxval(xmin) max_ymin = maxval(ymin) min_xmin = minval(xmin) min_ymin = minval(ymin) ! now update the pressure and density in the mhd code ! use the lfm grid to transfer the rcm information !-------------------------------------------------------------------- !FIXME: The following code converts pressure, density and mask from ! the RCM grid onto the LFM grid. This code should probably be moved ! to some "prepareExport" function, which should be called internally ! by the ExportMagnetosphere() subroutine. !-------------------------------------------------------------------- do k=1,nk_mid do j=1,nj_mid do i=1,ni_mid ctemp = 0.0 ptemp = 0.0 rhotemp = 0.0 mask = 0 bleed_PRESSURE_lfm(i,j,k) = ctemp bleed_RHO_lfm(i,j,k) = rhotemp bleed_MASK_lfm(i,j,k) = mask ! take each point on the midpoint of the 3-D LFM grid: xtrans = xmid_lfm(i,j,k) ytrans = ymid_lfm(i,j,k) ztrans = zmid_lfm(i,j,k) rtrans = sqrt(xtrans**2+ytrans**2+ztrans**2) if(rtrans < rtrans_min)CYCLE IF(use_ionosphere_interpolation)THEN call get_press_rho_from_ionos(idim,jdim,jwrap,& colat,aloct,imin_j,pressrcm, & densrcm,xtrans,ytrans,ztrans, & ptemp,rhotemp,mask,ifound,ierr2) ELSE call get_press_rho_from_equator(xtrans,ytrans,ztrans,& imin_grid,i_offset,ptemp,rhotemp,mask,ifound) END IF if(mask.eq.1.and.ptemp.eq.0.and.rhotemp.eq.0)then write(6,*)' mask =',mask write(6,*)' ptemp =',ptemp write(6,*)' rhotemp=',rhotemp write(6,*)' i j k =',i,j,k write(6,*)' ifound =',ifound end if ! all transfers are in SI units ! write transfer grid information and arrays when l_write_rcmu flag is ! on ! if(l_write_rcmu)call write_rcm_lfm_transfer_grid if(ifound.eq.2.and.mask.eq.1)then bleed_PRESSURE_lfm(i,j,k) = ptemp ! in Pa bleed_RHO_lfm(i,j,k) = rhotemp ! in kg/m^3 bleed_MASK_lfm(i,j,k) = mask end if end do end do end do ! if the locations are within 1 grid point of the boundary, then set the mask to zero if(avoid_boundaries)then bleed_MASK_lfm_temp = bleed_MASK_lfm do k=1,nk_mid do j=1,nj_mid do i=1,ni_mid im = max(1,i-1); ip = min(ni_mid,i+1) jm = max(1,j-1); jp = min(nj_mid,j+1) km = max(1,k-1); kp = min(nk_mid,k+1) if(k==1)km = nk_mid; if(k==nk_mid)kp = 1 istart = min(im,ip); iend = max(im,ip) jstart = min(jm,jp); jend = max(jm,jp) kstart = min(km,kp); kend = max(km,kp) bleed_MASK_lfm_temp(i,j,k) = minval(bleed_MASK_lfm(istart:iend,jstart:jend,kstart:kend)) end do end do end do ! reset the bleed mask to the new value bleed_MASK_lfm = bleed_MASK_lfm_temp end if !-------------------------------------------------------------------- !FIXME: The preceeding code should be moved inside ! ExportMagnetosphere... See previous "FIXME" comment for more info. !-------------------------------------------------------------------- ! Export magnetosphere variables in intermediate_grid module to LFM MHD. call ExportMagnetosphere() RETURN END SUBROUTINE tomhd ! !--------------------------------------- FUNCTION Rinter (x1,x2,x3,x4,dx,dy) USE Rcm_mod_subs, ONLY : rprec IMPLICIT NONE REAL(rprec), INTENT (IN) :: x1,x2,x3,x4,dx,dy REAL(rprec) :: Rinter ! ! x3-----------x4 ! | | | ! | | | ! +dy-+--------+ ! | | | ! | dx | ! | | | ! x1--+-------x2 ! ! Rinter = x4*dx*dy +x2*dx*(1.-dy)+ & x3*(1.-dx)*dy+x1*(1.-dx)*(1.-dy) ! RETURN END FUNCTION Rinter ! ! !***************************************************** ! l o c a t 3 !***************************************************** subroutine locat3(rval,pval,xval,yval,latdim,ltdim,bndi,r,p, & x,y,bi,bj,iout) ! purpose: find grid location of physical location given by ! (rval,pval) in cylindrical coordinates and (xval,yval) ! in gsm coordinates. ! ! version: 2.00 date: 6/21/94 USE Rcm_mod_subs, ONLY : iprec,rprec implicit none integer(iprec) :: latdim,ltdim,imax,juse real(rprec) :: r(latdim,ltdim),p(latdim,ltdim),x(latdim,ltdim) real(rprec) :: y(latdim,ltdim),bndi(ltdim) real(rprec) :: ycross(10) real(rprec) :: pi,dpmin,rval,pval,xval,yval,bi,bj real(rprec) :: p4,p5,dp real(rprec) :: x1,y1,x2,y2,x3,y3,x4,y4 real(rprec) :: a1,b1,c1,a2,b2,c2 real(rprec) :: u1,t1,tplus,tminus,u,uminus,t real(rprec) :: distmn,distsq integer(iprec) :: iout,j,ipt,jpt,i integer(iprec) :: ii,jj,jx,iuse,ifind pi=atan2(0.,-1.) imax=latdim ! find j value that has closest p value at i=imax dpmin=1.e6 juse=0 do 10 j=2,ltdim-2 ! make sure that pval and p(imax,j) are in same modulus call pfix(pval,p(imax,j),p(imax,j),p4,p5) dp=abs(pval-p4) if(dp.lt.dpmin) then juse=j dpmin=dp end if 10 continue ! write(6,*) 'closest j at imax=',juse if(juse.eq.0) then write(6,*) 'trouble in locat2 juse=',juse write(6,*)' x,y = ',xval,yval iout=9 return end if ! move out in i until r(i,juse).gt.rval while adjusting juse to ! continue to be closest j value to pval if(rval.le.r(imax,juse)) then iout=9 bi=999. bj=999. return end if iuse=imax 20 continue iuse=iuse-1 if(iuse.lt.2) then iout=2 return end if dpmin=1.e6 jx=0 do 25 j=juse-1,juse+1 jj=j if(j.lt.2) jj=j+ltdim-3 if(j.gt.ltdim-2) jj=j-ltdim+3 call pfix(pval,p(iuse,j),p(iuse,j),p4,p5) dp=abs(pval-p4) if(dp.lt.dpmin) then jx=jj dpmin=dp end if 25 continue if(jx.eq.0) then write(6,*) 'trouble in locat2 jx=',jx write(6,*)' x,y = ',xval,yval iout=9 return end if juse=jx ! write(6,*) 'iuse,juse',iuse,juse ! juse is closest j value to pval for i=iuse ! check if we are outside of bndy if(iuse.le.latdim-2) then if(float(iuse+2).lt.bndi(juse)) then iout=9 bi=999. bj=999. return end if end if ! check if r(iuse,juse) is still .lt. rval if(r(iuse,juse).lt.rval) go to 20 ! to move to next i ! look at nearby pts and determine closest grid pt distmn=1.e6 ipt=0 jpt=0 do 30 i=iuse,iuse+1 do 40 j=juse-1,juse jj=j if(j.lt.2) jj=j+ltdim-3 if(j.gt.ltdim-2) jj=j-ltdim+3 distsq=(x(i,jj)-xval)**2+(y(i,jj)-yval)**2 if(distsq.lt.distmn) then distmn=distsq ipt=i jpt=jj end if 40 continue 30 continue ! (ipt,jpt) is closest grid pt to (xval,yval) ! write(6,*) 'closest grid point ipt,jpt',ipt,jpt if(ipt.eq.0) then iout=4 return end if ! find grid cell containing (xval,yval)( ! first check 4 grid cells having closest point as a vertex do 100 i=ipt-1,ipt ii=max(i,1) ! added frt 2/98 ii=min(ii,imax-1) do 110 j=jpt-1,jpt jj=j if(jj.lt.1) jj=jj+ltdim-3 if(jj.gt.ltdim) jj=jj-ltdim+3 ! (ii,jj) gives location of lower left corner of cell to ! be checked call fndbox(ii,jj,xval,yval,x,y,latdim,ltdim,ifind) if(ifind.eq.1) then iuse=ii juse=jj go to 200 end if 110 continue 100 continue ! check more distant grid squares do 120 i=ipt-2,ipt+1 ii=max(i,1) ! added frt 2/98 ii=min(ii,imax-1) do 130 j=jpt-2,jpt+1 jj=j if(jj.lt.1) jj=jj+ltdim-3 if(jj.gt.ltdim) jj=jj-ltdim+3 if(ii.eq.ipt.and.j.eq.jpt-1) go to 130 if(ii.eq.ipt-1.and.j.eq.jpt-1) go to 130 if(ii.eq.ipt.and.j.eq.jpt) go to 130 if(ii.eq.ipt-1.and.j.eq.jpt) go to 130 call fndbox(ii,jj,xval,yval,x,y,latdim,ltdim,ifind) if(ifind.eq.1) then iuse=ii juse=jj go to 200 end if 130 continue 120 continue ! unable to find proper grid square. set bi,bj equal to closest ! grid point. set iout=-1 and return. if(float(ipt).lt.bndi(jpt)) then iout=9 bi=999. bj=999. return else iout=-1 bi=float(ipt) bj=float(jpt) return end if 200 continue ! proper grid square for interpolation is defined by (iuse,juse), ! (iuse,juse+1),(iuse+1,juse+1), and (iuse+1,juse). ! write(6,*) 'xval,yval,iuse,juse',xval,yval,iuse,juse if((float(iuse).lt.bndi(juse)).or.(float(iuse).lt.bndi(juse+1)))& then iout=9 bi=999. bj=999. return end if iout=0 x1=x(iuse,juse) x2=x(iuse,juse+1) x3=x(iuse+1,juse+1) x4=x(iuse+1,juse) y1=y(iuse,juse) y2=y(iuse,juse+1) y3=y(iuse+1,juse+1) y4=y(iuse+1,juse) ! write(6,*) 'x1,x2,x3,x4,y1,y2,y3,y4',x1,x2,x3,x4,y1,y2,y3,y4 a1=(x4-x1)*(y2-y3)+(x2-x3)*(y1-y4) b1=(xval-x1)*(y1-y2+y3-y4)+(x1-x4)*(y2-y1)- & (yval-y1)*(x1-x2+x3-x4)-(y1-y4)*(x2-x1) c1=(xval-x1)*(y2-y1)-(yval-y1)*(x2-x1) a2=(x2-x1)*(y4-y3)+(x4-x3)*(y1-y2) b2=(xval-x1)*(y1-y2+y3-y4)+(x1-x2)*(y4-y1)- & (yval-y1)*(x1-x2+x3-x4)-(y1-y2)*(x4-x1) c2=(xval-x1)*(y4-y1)-(yval-y1)*(x4-x1) ! write(6,*) 'a1,b1,c1',a1,b1,c1 ! write(6,*) 'a2,b2,c2',a2,b2,c2 u1=-c1/b1 call quadf(a1,b1,c1,u,uminus) ! ux=u1-(a1*u1**2+b1*u1+c1)/(2.*a1*u1+b1) t1=-c2/b2 call quadf(a2,b2,c2,tplus,t) ! tx=t1-(a2*t1**2+b2*t1+c2)/(2.*a2*t1+b2) ! write(6,*) 'u1,t1,ux,tx',u1,t1,ux,tx bi=float(iuse)+u bj=float(juse)+t ! 500 continue ! RETURN END SUBROUTINE Locat3 ! ! !********************************** pfix subroutine pfix(p1,p2,p3,p4,p5) ! copyright rice university, 1993 ! purpose: subroutine to adjust modulus of p2 and p3 to match p1. ! on output p4 corresponds to p2, p5 to p3 ! ! version 1.00 date: 09.04.89 ! 2.00 msm delivery version 02.02.93 ! 2.10 debug of p's near zero 08.09.93 ! ! programmer: r.w. spiro USE Rcm_mod_subs, ONLY : rprec implicit none real(rprec) :: p1,p2,p3,p4,p5,pi pi=atan2(0.,-1.) ! filter for msm grid locations with p's near zero if (abs(p2).lt.(1.e-5)) p2 = 0. if (abs(p3).lt.(1.e-5)) p3 = 0. p4=p2 p5=p3 if(abs(p2-p3).gt.pi) then if(p2.lt.p3) then p4=p2+2.*pi else p5=p3+2.*pi end if if(p1.lt.pi) then p4=p4-2.*pi p5=p5-2.*pi end if end if if(abs(p2-p3).gt.pi) then if(p2.lt.p3) then p4=p2+2.*pi else p5=p3+2.*pi end if if(p1.lt.pi) then p4=p4-2.*pi p5=p5-2.*pi end if end if return end !******************************************** ! q u a d f !******************************************** subroutine quadf(a,b,c,root1,root2) ! roots of a quadratic eqn USE Rcm_mod_subs, ONLY : rprec implicit none real(rprec) :: a,b,c,root1,root2 if(abs(a).gt.1.e-16) then root1=(-b+sqrt(b**2-4.*a*c))/(2.*a) root2=(-b-sqrt(b**2-4.*a*c))/(2.*a) else root1=-c/b root2=root1 end if return end !*********************************************************** ! x i n g !*********************************************************** subroutine xing(x1,x2,y1,y2,xval,yxx,ierr) USE Rcm_mod_subs, ONLY : iprec,rprec implicit none real(rprec) :: x1,x2,y1,y2,xval,yxx integer(iprec) :: ierr ! purpose: subroutine to determine crossing point of line from ! (x1,y1,) to (x2,y2) with line x=xval. ! if the two lines are parallel (or nearly so) ierr is set ! equal to 1 and the subroutine returns; otherwise, yxx ! is the crossing point and ierr is set to zero. ! date: 07.15.94 if(abs(x1-x2).lt.1.e-4) then ierr=1 else yxx=y1+(xval-x1)*(y2-y1)/(x2-x1) ierr=0 end if return end !********************************************************* ! f n d b o x !********************************************************* subroutine fndbox(ii,jj,xval,yval,x,y,latdim,ltdim, & ifind) ! version 1.00 07.18.94 ! purpose: subroutine to determine if grid square with (ii,jj) ! as lower left corner contains point (xval,yval). ! if so, set ifind=1; otherwise, ifind=0 USE Rcm_mod_subs, ONLY : iprec,rprec implicit none integer(iprec) :: ii,jj,latdim,ltdim,ifind real(rprec) :: xval,yval real(rprec) :: xmax,xmin,ymax,ymin real(rprec) :: x1,y1,x2,y2 real(rprec) :: temp,yxx real(rprec) :: x(latdim,ltdim),y(latdim,ltdim) real(rprec) :: ycross(10) integer(iprec) :: ix,jx,ncross,i1,i2,j1,j2,n,nx,iside,ierr ifind=0 ! find extremes of grid cell xmax=-1.e6 xmin=1.e6 ymax=-1.e6 ymin=1.e6 ! do 120 ix=ii,ii+1 do 130 jx=jj,jj+1 if(x(ix,jx).gt.xmax) xmax=x(ix,jx) if(x(ix,jx).lt.xmin) xmin=x(ix,jx) if(y(ix,jx).gt.ymax) ymax=y(ix,jx) if(y(ix,jx).lt.ymin) ymin=y(ix,jx) 130 continue 120 continue if(xval.lt.xmin.or.xval.gt.xmax.or.yval.lt.ymin.or. & yval.gt.ymax) then ifind=0 return end if ! otherwise, find intersection of x=xval line with each ! side of grid cell if intersections are inside of extremes. ncross=0 do 140 iside=1,4 if(iside.eq.1) then i1=ii i2=ii j1=jj j2=jj+1 else if(iside.eq.2) then i1=ii i2=ii+1 j1=jj+1 j2=jj+1 else if(iside.eq.3) then i1=ii+1 i2=ii+1 j1=jj+1 j2=jj else i1=ii+1 i2=ii j1=jj j2=jj end if ! find value of y where line x=xval crosses side. call this ! value yxx. ierr.ne.0 means x=xval parallels side x1=x(i1,j1) y1=y(i1,j1) x2=x(i2,j2) y2=y(i2,j2) call xing(x1,x2,y1,y2,xval,yxx,ierr) if(ierr.ne.0) go to 140 ! check if yxx is between y1 and y2 if((yxx.le.y1.and.yxx.ge.y2).or. & (yxx.le.y2.and.yxx.ge.y1)) then ! write(6,*) '** ii,jj,iside',ii,jj,iside, '*****' ! write(6,*) 'x1,y1,x2,y2',x1,y1,x2,y2 ! write(6,*) 'xval,yval,yxx,ierr',xval,yval,yxx,ierr ! valid crossing; count it ncross=ncross+1 ycross(ncross)=yxx ! write(6,*) 'ncross,ycross(ncross)',ncross,yxx end if 140 continue if(ncross.ne.2.and.ncross.ne.4) then ifind=0 return else ! order values in ycross using piksrt algorithm per ! press et al., numerical recipes, p.227. do 150 n=2,ncross temp=ycross(n) do 152 nx=n-1,1,-1 if(ycross(nx).le.temp) go to 154 ycross(nx+1)=ycross(nx) 152 continue nx=0 154 ycross(nx+1)=temp 150 continue if(yval.ge.ycross(1).and.yval.le.ycross(2)) then ! (xval,yval) lays within the grid square defined ! by (ii,jj),(ii,jj+1),(ii+1,jj),(ii+1,jj) ifind=1 return else if(ncross.eq.4.and. & (yval.ge.ycross(3).and.yval.le.ycross(4))) then ifind=1 return else ifind=0 return end if end if ! return end !----------------------------------------------------------- subroutine writeu_2d(idim,jdim,imax,jmax,file,array,x,y,rec) ! routine to writeout array unformatted at some interval ! 8/08? frt USE Rcm_mod_subs, ONLY : iprec,rprec implicit none integer(iprec) :: idim,jdim integer(iprec) :: imax,jmax real(rprec) :: array(idim,jdim) real(rprec) :: x(idim),y(jdim) integer(iprec) :: rec,i character (LEN=*) :: file ! if(rec.eq.1)then if(rec.le.2)then open(80,file=file,status='unknown',form='unformatted') write(80)idim,jdim,imax,jmax write(80)rec,array,x,y close(80) else open(80,file=file,status='unknown',form='unformatted' & ,position='append') write(80)rec,array,x,y close(80) end if return end !--------------------------------------------pressdens ! pressure function initial setup for setupd.f ! follows a field line until it gets to the equatorial ! plane then assigns a pressure according to the empirical result ! work in gsm ! modified 9/26/03 to also return density - frt ! added ifound flag to see if you cross the equatorial plane subroutine find_eqcross(x,y,z,x_new,y_new,z_new,ifound,ierr) USE intermediate_grid, ONLY : xmin_ig,xmax_ig,& ymin_ig,ymax_ig,& zmax_ig,zmin_ig USE Rcm_mod_subs, ONLY : iprec,rprec implicit none real(rprec) :: x,y,z real(rprec) :: x_old,y_old,z_old,v_old real(rprec) :: x_new,y_new,z_new,v_new real(rprec) :: b real(rprec) :: er,dir,h real(rprec) :: xmax,xmin real(rprec) :: ymax,ymin real(rprec) :: zmax,zmin real(rprec) :: r2min real(rprec) :: dx,dy,dz,bt,bt0 integer(iprec) :: i,ierr,ifound real(rprec), parameter :: hmax = 1.0 real(rprec), parameter :: hmin = 0.01 real(rprec), parameter :: rmin = 1.0 integer(iprec), parameter :: imax = 5000 ! changed 10/19/2004, SS data rmin /3.5/ x_old = x y_old = y z_old = z x_new=x y_new=y z_new=z r2min = rmin**2 h = 1.0 dir = -1.0 if(z.lt.0.0)dir=1.0 er = 0.005 ifound = 0 do i=1,imax ! if((x_old**2+y_old**2+z_old**2).lt.r2min)go to 1 call tracf2(x_old,y_old,z_old,v_old,& x_new,y_new,z_new,v_new,er,dir,h,ierr) if(ierr.lt.0)return ! write(100,*)i,x_new,y_new,z_new if(h.ge.hmax)h=hmax if(h.le.hmin)h=hmin if(x_new.gt.xmax_ig.or.x_new.lt.xmin_ig)return if(y_new.gt.ymax_ig.or.y_new.lt.ymin_ig)return if(z_new.gt.zmax_ig.or.z_new.lt.zmin_ig)return if((x_new**2+y_new**2+z_new**2).lt.r2min)return if(z_new*z_old.le.0.0)then ! reset to z=0 if(z_new-z_old.ne.0)then x_new = (x_new*z_old - x_old*z_new)/(z_old-z_new) y_new = (y_new*z_old - y_old*z_new)/(z_old-z_new) z_new = 0.0 end if ifound = 1 return end if x_old = x_new y_old = y_new z_old = z_new end do write(6,*)' limit in i reached in pressure, aborting' write(6,*)' first point at',x,y,z write(6,*)' last point at',x_new,y_new,z_new return end subroutine find_eqcross !-----------------------inter subroutine inter2(nx,ny,nxp,nyp,xp,yp, & press,ipflag,x,y,p,ierr) USE Rcm_mod_subs, ONLY : iprec,rprec implicit none integer(iprec) :: nx,ny,nxp,nyp integer(iprec) :: ierr real(rprec) :: xp(nx),yp(ny),press(nx,ny),x,y,p integer(iprec) :: ipflag(nx,ny) integer(iprec) :: i,j,ig,jg real(rprec) :: dx,dy real(rprec) :: a1,a2,a3,a4 real(rprec) :: nonz_ave p = 0.0 do i=1,nxp-1 if(x.ge.xp(i).and.x.lt.xp(i+1))then ig = i dx = x - xp(i) go to 1 end if end do return 1 do j=1,nyp-1 if(y.ge.yp(j).and.y.lt.yp(j+1))then jg = j dy = y - yp(j) go to 2 end if end do return 2 continue ! i,j+1---+-----i+1,j+1 ! | | | ! | a3 | a4 | ! +-----+--------+ ! | | | ! | a1 dy a2 | ! | | | ! +-dx--+------i+1,j a1 = dx*dy a2 = (xp(ig+1)-x)*dy a3 = dx*(yp(jg+1)-y) a4 = (xp(ig+1)-x)*(yp(jg+1)-y) if(a1.lt.0.0.or.a2.lt.0.0.or.a3.lt.0.0.or.a4.lt.0.0)then write(6,*)' error in inter ' write(6,*)' a1,a2,a3,a4 =',a1,a2,a3,a4 ierr = -2 return end if ! if any of the pressure terms are zero then zero a's if(ipflag(ig,jg).eq.0)a4 = 0.0 if(ipflag(ig+1,jg).eq.0)a3=0.0 if(ipflag(ig,jg+1).eq.0)a2=0.0 if(ipflag(ig+1,jg+1).eq.0)a1=0.0 p = press(ig,jg)*a4 + & press(ig+1,jg)*a3 + & press(ig,jg+1)*a2 + & press(ig+1,jg+1)*a1 if(a1+a2+a3+a4.gt.0.0)then p = p/(a1+a2+a3+a4) else p = nonz_ave(press(ig,jg)*ipflag(ig,jg), & press(ig+1,jg)*ipflag(ig+1,jg), & press(ig,jg+1)*ipflag(ig,jg+1), & press(ig+1,jg+1)*ipflag(ig+1,jg+1)) end if ! if(p.le.0.0)then ! write(6,*)' !!!!!!!!!!!!!!debug point in pressure p <=0' ! end if return end !********************************************* real function nonz_ave(p1,p2,p3,p4) implicit none ! takes the average of 4 non-zero terms ! 4/96 frt real p1,p2,p3,p4 real sum,tot sum = 0.0 tot = 0.0 nonz_ave = 0.0 if(p1.ne.0.0)then sum=sum+p1 tot = tot + 1.0 end if if(p2.ne.0.0)then sum=sum+p2 tot = tot + 1.0 end if if(p3.ne.0.0)then sum=sum+p3 tot = tot + 1.0 end if if(p4.ne.0.0)then sum=sum+p4 tot = tot + 1.0 end if if(tot.ne.0.0)nonz_ave=sum/tot return end !--------------------------------------------get_press_rho_from_ionos ! pressure function initial setup for setupd.f ! follows a field line until it gets to the ionosphere ! written 2/98 frt ! p(x) = 5.3e-8 exp(-x/2.12) x is in Re ! work in gsm subroutine get_press_rho_from_ionos(idim,jdim,jwrap,colat,aloct,imin_j,pressrcm, & densrcm,x,y,z,p,rho,mask,ifound,ierr) !---------------------------------------------------------------------- ! inputs: ! idim, jdim - dimensions of rcm grid ! colat and aloct - colatitude and long in radians ! press - pressure distibution on rcm grid ! imin_j - first point in rcm grid that is inside region ! x,y,z - location on physical space for which pressure is ! to be found ! p - returned pressure ! ierr - error flag ! -1 means abort program due to a serious error ! 0 means a viable pressure was found ! 1 means a field line not connected to ionosphere ! written 2/19/98 frt !---------------------------------------------------------------------- USE intermediate_grid, ONLY : xmin_ig,xmax_ig,& ymin_ig,ymax_ig,& zmax_ig,zmin_ig USE Rcm_mod_subs, ONLY : iprec,rprec USE constants, ONLY : mass_proton implicit none integer(iprec) ::idim,jdim,jwrap real(rprec) :: colat(idim),aloct(idim) real(rprec) :: pressrcm(idim,jdim),rhorcm(idim,jdim) real(rprec) :: densrcm(idim,jdim) integer(iprec) :: imin_j(jdim) real(rprec) :: x,y,z real(rprec) :: x_old,y_old,z_old,v_old real(rprec) :: x_new,y_new,z_new,r_new,v_new real(rprec) :: rr real(rprec) :: xout,yout,zout real(rprec) :: p,b real(rprec) :: er,er0,dir,h real(rprec) :: hmax,hmin real(rprec) :: rmin,r2min real(rprec) :: rho real(rprec) :: dx,dy,dz,bt,bt0 real(rprec) :: p1,p2,p3 real(rprec), parameter :: grr = 1.0 ! radius of earth integer(iprec) :: i,imax,mask integer(iprec) :: ierr integer(iprec),intent(out) :: ifound data hmax /0.1/ data hmin /0.01/ data imax /4000/ data b /0.4738/ data rmin /2.0/ ierr = 0 ifound = 2 x_old = x y_old = y z_old = z v_old = 0.0 r2min = rmin**2 p = 0.0 rho = 0.0 mask = 0 rhorcm = mass_proton * densrcm h = 1.0 dir = 1.0 er0 = 0.01 do i=1,imax rr = x_old**2+y_old**2+z_old**2 er = er0*(rr/100.) ! make the error smaller near the earth call tracf2(x_old,y_old,z_old,v_old,& x_new,y_new,z_new,v_new,er,dir,h,ierr) if(ierr.lt.0)return ! write(100,*)i,x_new,y_new,z_new h=amin1(h,hmax) h=amax1(h,hmin) if(x_new.gt.xmax_ig.or.x_new.lt.xmin_ig)go to 2 if(y_new.gt.ymax_ig.or.y_new.lt.ymin_ig)go to 2 if(z_new.gt.zmax_ig.or.z_new.lt.zmin_ig)go to 2 r_new = sqrt(x_new**2+y_new**2+z_new**2) ! hit the ionosphere if(r_new.le.rmin)go to 1 x_old = x_new y_old = y_new z_old = z_new v_old = v_new end do write(6,*)' limit in i reached in pressure, aborting, i =',i write(6,*)' first point at',x,y,z write(6,*)' last point at + v',x_new,y_new,z_new,v_new ierr = -1 return 1 rho = r_new ! interpolate to r=grr call trans_dipt(0.0,x_new,y_new,z_new,grr,xout,yout,zout) call interi(idim,jdim,jwrap,colat,aloct,imin_j, & pressrcm,xout,yout,zout,p,ierr) call interi(idim,jdim,jwrap,colat,aloct,imin_j, & rhorcm,xout,yout,zout,rho,ierr) mask = 1 ifound = 2 if(ierr == 1)then ! outside the modeling region p = 0.0 rho = 0.0 mask = 0 ifound = 1 return endif if(ierr.lt.0)return ! if all else fails use the old value of pressure ! if(p.le.0.0)then ! call bgrid(2,x,y,z,p1,p2,p3,ierr) ! if(ierr.lt.0)return ! p=p1 ! end if return 2 ierr = 1 mask = 0 ifound = 0 return end subroutine get_press_rho_from_ionos !-----------------------interi---------------------------------------- subroutine interi(idim,jdim,jwrap,colat,aloct, & min_j,pressi,x,y,z,p,ierr) !--------------------------------------------------------------------- ! this subroutine interpolates in the ionosphere ! writtend 2/20/98 frt ! inputs: ! idim, jdim - dimensions of rcm grid ! colat and aloct - colatitude and long in radians ! min_j - loaction of the first grid point inside bndy ! press - pressure distibution on rcm grid ! x,y,z - location on physical space for which pressure is ! to be found ! p - returned pressure ! ierr - error flag ! -1 means abort program due to a serious error !--------------------------------------------------------------------- USE Rcm_mod_subs, ONLY : iprec,rprec implicit none integer(iprec) :: idim,jdim,jwrap integer(iprec) :: ierr real(rprec) :: colat(idim,jdim),aloct(idim,jdim) real(rprec) :: pressi(idim,jdim),x,y,z,p integer(iprec) :: min_j(jdim) integer(iprec) :: i,j,ig,jg real(rprec) :: colatp,aloctp,r real(rprec) :: dlat,dlong real(rprec) :: a1,a2,a3,a4 real(rprec) :: p1,p2,p3 real(rprec) :: nonz_ave real(rprec) :: pi p = 0.0 ierr = 0 pi = acos(-1.0) ! find colatp and aloctp r=sqrt(x**2+y**2+z**2) if(r.le.0.0)then write(6,*)' r = 0 in interi' ierr = -1 return end if if(z.lt.0.0)then write(6,*)' error: z lt 0 in interi' ierr = -1 return end if colatp = acos(z/r) if(y.eq.0.0.and.x.eq.0.0)then aloctp = 0.0 else aloctp = atan2(y,x) end if if(aloctp.lt.0.0)aloctp = aloctp + 2.0*pi ! find grid indices if(colatp.lt.colat(1,1).or.colatp.gt.colat(idim,1))then ierr = 1 return end if do i=1,idim-1 if(colatp.ge.colat(i,1).and.colatp.lt.colat(i+1,1))then ig = i dlat = colatp - colat(i,1) go to 1 end if end do ierr = -1 write(6,*)' ig not found in interi' write(6,*)' colatp,aloctp =',colatp/pi*180.,aloctp/pi*180. return 1 do j=1,jdim-2 if(aloctp.ge.aloct(1,j).and.aloctp.lt.aloct(1,j+1))then jg = j dlong = aloctp - aloct(1,j) go to 2 end if end do ! last resort if the point is between aloct(1,jdim-1) and 2pi if(aloctp.ge.aloct(1,jdim-1).and.aloctp.le.2.0*pi)then jg = jdim - 1 dlong = aloctp - aloct(i,jdim-1) go to 2 end if write(6,*)' jg not found in interi' write(6,*)' colatp,aloctp =',colatp/pi*180.,aloctp/pi*180. ierr = -1 return 2 continue ! check to see if we are outside the modeling region if(ig.lt.min_j(jg))then ierr = 1 return end if ! i,j+1---+-----i+1,j+1 ! | | | ! | a3 | a4 | ! +-----+--------+ ! | | | ! | a1 dy a2 | ! | | | ! +-dx--+------i+1,j if(jg.lt.jdim-1)then a1 = dlat*dlong a2 = (colat(ig+1,1)-colatp)*dlong a3 = dlat*(aloct(1,jg+1)-aloctp) a4 = (colat(ig+1,1)-colatp)*(aloct(1,jg+1)-aloctp) else a1 = dlat*dlong a2 = (colat(ig+1,1)-colatp)*dlong a3 = dlat*(2.0*pi-aloctp) a4 = (colat(ig+1,1)-colatp)*(2.0*pi-aloctp) end if if(a1.lt.0.0.or.a2.lt.0.0.or.a3.lt.0.0.or.a4.lt.0.0)then write(6,*)' error in inter ' write(6,*)' a1,a2,a3,a4 =',a1,a2,a3,a4 ierr = -1 return end if ! if any of the pressure terms are zero then zero a's if(pressi(ig,jg).le.0.0) a4=0.0 if(pressi(ig+1,jg).le.0.0) a3=0.0 if(pressi(ig,jg+1).le.0.0) a2=0.0 if(pressi(ig+1,jg+1).le.0.0)a1=0.0 p = pressi(ig,jg)*a4 + & pressi(ig+1,jg)*a3 + & pressi(ig,jg+1)*a2 + & pressi(ig+1,jg+1)*a1 if(a1+a2+a3+a4.gt.0.0)then p = p/(a1+a2+a3+a4) else p = nonz_ave(pressi(ig,jg), pressi(ig+1,jg), & pressi(ig,jg+1),pressi(ig+1,jg+1)) end if if(p.le.0.0)then ! call bgrid(2,x,y,z,p1,p2,p3) ! p=p1 ! write(6,*)' !!!!!!!!!!!!!!debug point in pressure p <=0' end if return end !--------------------------------------------------- subroutine get_press_rho_from_equator(xtrans,ytrans,ztrans,& imin_grid,i_offset,ptemp,rhotemp,mask,ifound) ! gets the pressure and rho from the rcm equatorial mapped grid ! ifound =0 - not finding a mapping ! ifound =1 - found a mapping, but outside the rcm's region ! ifound =2 - found a valid result USE constants, ONLY : mass_proton USE RCM_mod_subs, ONLY: iprec,rprec,isize,jsize,jwrap,& pmin,rmin,pressrcm,densrcm,bndloc,xmin,ymin,vm IMPLICIT NONE REAL(rprec), intent(in) :: xtrans,ytrans,ztrans REAL(rprec), intent(out) :: rhotemp,ptemp INTEGER(iprec), intent(out) :: mask REAL(rprec) :: xn,yn,zn,rval,pval REAL(rprec) :: max_xmin,max_ymin,min_xmin,min_ymin REAL(rprec) :: bbi,bbj,di,dj REAL(rprec) :: pi REAL(rprec) :: rinter REAL(rprec) :: rhorcm(isize,jsize) INTEGER(iprec), intent(in) :: imin_grid,i_offset INTEGER(iprec), intent(out) :: ifound INTEGER(iprec) :: iout,ig,jg,idim,jdim,ierr2 rhorcm = densrcm * mass_proton pi = acos(-1.0) max_xmin = maxval(xmin) max_ymin = maxval(ymin) min_xmin = minval(xmin) min_ymin = minval(ymin) idim = isize jdim = jsize ifound = 0 ! means not found ! and trace from this point to the equtorial plane: CALL Find_eqcross (xtrans,ytrans,ztrans,xn,yn,zn,ifound,ierr2) ! If could not find an equatorial plane crossing (open field line?), ! then forget about this point and move on to the next grid point: IF (ifound == 0) return ! Otherwise, have equatorial crossing at (xn,yn,zn). ! If this crossing is outside the RCM modeling region, forget about ! this point and move on to the next grid point: IF (xn >= max_xmin .or. & xn <= min_xmin .or. & yn >= max_ymin .or. & yn <= min_ymin) RETURN ! Otherwise, compute the radius R and local-time angle P: rval = SQRT (xn**2+yn**2) IF (xn == 0.0 .AND. yn == 0.0) THEN pval = 0.0 ELSE pval = ATAN2 (yn,xn) END IF IF (pval < 0.0) pval = pval+2.0*pi ! Find non-integer indices for location of this point on the RCM grid: CALL Locat3 (rval, pval, xn, yn, idim, jdim, & bndloc, rmin, pmin, xmin, ymin, bbi, bbj, iout) IF (iout /= 0) return IF (bbi < 1.0 .OR. bbj < 1.0) return ig = INT (bbi) jg = INT (bbj) IF (ig <= imin_grid+ i_offset) return IF (ig > idim) return ! distances: di = bbi - REAL(ig) dj = bbj - REAL(jg) ! If any of the four corners of the RCM grid cell containing the eq.plane ! crossing is on a field line that is marked open, don't do update, ! and move on. It is technically possible to have (xmin,ymin) values ! in the RCM for open field lines, but not via field-line tracing. IF (vm(ig,jg) <= 0.0) return IF (vm(ig,jg+1) <= 0.0) return IF (vm(ig+1,jg) <= 0.0) return IF (vm(ig+1,jg+1) <= 0.0) return ! Finally, interpolate pressure and density for the location: ptemp = Rinter (pressrcm(ig,jg), pressrcm(ig+1,jg), & pressrcm(ig,jg+1), pressrcm(ig+1,jg+1),di,dj) rhotemp = Rinter (rhorcm(ig,jg), rhorcm(ig+1,jg), & rhorcm(ig,jg+1), rhorcm(ig+1,jg+1),di,dj) mask = 1 ifound =2 RETURN end subroutine get_press_rho_from_equator !--------------------------------------------------- ! subroutine to transform to r=r from a position in ! space a xin,yin,zin to xout you zout ! includes dipole tilt ! 2/97 frt !--------------------------------------------------- subroutine trans_dipt(tilt,xin,yin,zin,r,xout,yout,zout) USE Rcm_mod_subs, ONLY : rprec USE rice_housekeeping_module implicit none real(rprec) :: xin,yin,zin real(rprec) :: xrot,yrot,zrot real(rprec) :: xout,yout,zout real(rprec) :: r real(rprec) :: r0,phi0,theta0 real(rprec) :: phi,theta real(rprec) :: factor real(rprec) :: tilt real(rprec) :: dpi ! compute r0,phi0,theta0 xout = 0.0 yout = 0.0 zout = 0.0 dpi = acos(-1.00)/180.0 xrot = xin*cos(tilt*dpi)-zin*sin(tilt*dpi) yrot = yin zrot = zin*cos(tilt*dpi)+xin*sin(tilt*dpi) ! write(6,*)' xrot,yrot,zrot =',xrot,yrot,zrot if(xrot.eq.0.0.and.yrot.eq.0.0)then xout = 0.0 yout = 0.0 zout = r*sign(1.0,zin) return end if r0 = sqrt(xrot**2+yrot**2+zrot**2) if(r0.ne.0.0)then theta0 =acos(zrot/r0) phi0 = atan2(yrot,xrot) ! write(6,*)' theta0=',theta0/3.14*180. ! write(6,*)' phi0=',phi0/3.14*180. else return end if ! now transform factor = sin(theta0)*sqrt(r/r0) ! no solution if(abs(factor).gt.1.0)return theta = asin(factor) xrot = r*cos(phi0)*sin(theta) yrot = r*sin(phi0)*sin(theta) zrot = r*cos(theta) xout = xrot*cos(tilt*dpi)+zrot*sin(tilt*dpi) yout = yrot zout = zrot*cos(tilt*dpi)-xrot*sin(tilt*dpi) ! if(L_write_vars_debug)then ! write(6,*)' xrot,yrot,zrot,r =',xrot,yrot,zrot, ! & sqrt(xrot**2+yrot**2+zrot**2) ! write(6,*)' xout,yout,zout,r =',xout,yout,zout, ! & sqrt(xout**2+yout**2+zout**2) ! end if return end !-----------------point----------------- subroutine point(x,y,z,bx,by,bz,bb) USE Rcm_mod_subs, ONLY: iprec,rprec implicit none real(rprec) :: x,y,z,bx,by,bz,bb real(rprec) :: xtrans,ytrans,ztrans,b1,b2,b3 integer(iprec) :: ierr xtrans = x ytrans = y ztrans = z call bfield_lfm(xtrans,ytrans,ztrans,b1,b2,b3,bb,ierr) if(bb.ne.0.0)then bx = b1 by = b2 bz = b3 else bx = 0.0 by = 0.0 bz = 0.0 end if return end !--------------------------------------------------- subroutine write_rcmu (record,offset) ! routine to write out rcm data in one unformated data file ! based on rcmu.dat 10/07 frt USE RCM_MOD_SUBS IMPLICIT NONE ! offset adds a time to itime to offset it INTEGER(iprec), INTENT(IN) :: record,offset INTEGER(iprec) :: itime IF (record == 1) THEN ! first time we write, record = 1 OPEN (LUN, FILE=rcmdir//'rcmu.dat', status = 'replace', & FORM = 'UNFORMATTED') ELSE OPEN (LUN, FILE=rcmdir//'rcmu.dat', status = 'old', & FORM = 'UNFORMATTED', POSITION='APPEND') END IF itime = label%intg(6) + offset write(6,*)'--->writing rcmu.dat at record =',record,' time =',itime WRITE (LUN) itime, isize, jsize, kcsize WRITE (LUN) alamc, etac, ikflavc WRITE (LUN) xmin,ymin,0*rmin, & sin(colat)*cos(aloct), & sin(colat)*sin(aloct),cos(colat), & vm,v_avg,eeta,birk_avg,bmin !,dsob3 close(LUN) return end subroutine write_rcmu subroutine carpenter(L,density) ! based on an approximation of carpenter, jgr 1992, figure 4 ! JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 97, NO. A2, PAGES 1097-1108, FEBRUARY 1, 1992 ! returns values in ples/cc USE Rcm_mod_subs, ONLY : rprec implicit none real(rprec),intent(in) :: L real(rprec),intent(out) :: density real(rprec), parameter :: nmax = 8022. ! ple/cc density = nmax / (max(1.0,L))**3.1 return end subroutine carpenter subroutine gallagher(L,density) ! approx based on gallagher et al, figure 1 ! JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 105, NO. A8, PAGES 18,819-18,833, AUGUST 1, 2000 ! returns values in ples/cc USE Rcm_mod_subs, ONLY : rprec implicit none real(rprec),intent(in) :: L real(rprec),intent(out) :: density real(rprec), parameter :: L0 = 4.5 real(rprec), parameter :: alpha = 10. real(rprec), parameter :: a1 = -0.25 real(rprec), parameter :: b1 = 2.4 real(rprec), parameter :: a2 = -0.5 real(rprec), parameter :: b2 = 4.5 real ::f,q f = 0.5*(1.0+tanh(alpha*(L-L0))) q = f*(a1*L + b1) + (1.0-f)*(a2*L + b2) density = 10.**q return end subroutine gallagher SUBROUTINE write_rcm_lfm_transfer_grid ! 7/09 frt ! writes in mks unis USE LFM_TRANSFERGRID USE CONSTANTS, ONLY : radius_earth_m USE Rcm_mod_subs, ONLY : iprec,rprec IMPLICIT none INTEGER(iprec) :: i,j,k open(unit=20,file='lfmgrid_data.dat',status='unknown',recl=120) write(20,*)ni_mid,nj_mid,nk_mid write(6,*)' lfmgrid_data , ni_mid =',ni_mid,' nj_mid=',& nj_mid,' nk_mid=',nk_mid ! use SI units to transfer do i=1,ni_mid do j=1,nj_mid do k=1,nk_mid write(20,*)xmid_lfm(i,j,k)*radius_earth_m,& ymid_lfm(i,j,k)*radius_earth_m,& zmid_lfm(i,j,k)*radius_earth_m,& bleed_rho_lfm(i,j,k),& bleed_pressure_lfm(i,j,k),& bleed_mask_lfm(i,j,k) end do end do end do close(20) return end subroutine write_rcm_lfm_transfer_grid