#include "dims.h" ! ! Routines in this file (inline.f) are inlined into their ! calling routines with the -Oinlinefrom=inline.f option ! to f90. Note that if there are any changes to this source, ! then the "calling" routines (i.e., the routines into which ! the changed routines are inlined) must also be recompiled ! for the changes to take effect. ! !----------------------------------------------------------------------- SUBROUTINE GRDINT(AM,AG,IG,JG,WT,LG,LM,IMAXG,JMAXG,IMAXM,JMAXM, | JM,iprint) implicit none C **** C **** TRANSFORM SCALAR FIELD GIVEN ON GEOGRAPHIC GRID TO C **** GEOMAGNETIC GRID USING INDICES AND WEIGHTS GENERATED C **** BY GRDSTM. C **** C **** INPUT: C **** AG(LG,1) IS 2-DIM SCALAR FIELD TO BE TRANSFORMED TO C **** GEOMAGNETIC GRID. C **** NB: PERIODIC POINT IS REPEATED SO THAT C **** AG(IMAXG+1,JG) = AG(1,JG) C **** (LM.GE.IMAXG+1) C **** IMAXG, JMAXG DIMENSIONS OF GEOGRAPHIC GRID C **** IMAXM, JMAXM DIMENSIONS OF GEOMAGNETIC GRID C **** LG IS FIRST DIMENSION OF AG IN CALLING PROGRAM C **** (LG.GE.IMAXG+1) C **** IG(LM,1) GIVES ROUNDED DOWN GEOGRAPHIC LONGITUDE INDEX C **** FOR EACH GEOMAGNETIC GRID POINT C **** JG(LM,1) GIVES ROUNDED DOWN GEOGRAPHIC LATITUDE INDEX C **** FOR EACH GEOMAGNETIC GRID POINT C **** WT(4,LM,1) WEIGHTS FOR 4 CORNERS OF EACH GEOMAGNETIC C **** GRID ELEMENT C **** LM IS FIRST DIMENSION OF ARRAYS AM, IG, JG, WT IN C **** CALLING PROGRAM C **** (LM.GE.IMAXM) C **** IEM IS FLAG INDICATING TYPE OF TRANSFORMATION C **** IEM = 1 FOR OUPUT FIELDS AT DIP EQUATOR C **** IEM = 0 FOR GLOBAL OUTPUT FIELDS C **** NOTE: IG, JG, WT ARE PRODUCED BY PREVIOUS CALL TO C **** GRDSTM. C **** C **** OUTPUT: C **** AM(LM,1) 2-DIM SCALAR FIELD TRANSFORMED TO GEOMAGNETIC C **** LATITUDE/LONGITUDE GRID C **** ! ! Args: integer,intent(in) :: lg,lm,imaxg,jmaxg,imaxm,jmaxm,jm,iprint real,intent(out) :: AM(LM,1) real,intent(in) :: AG(LG,1),WT(4,LM,1) integer,intent(in) :: IG(LM,1),JG(LM,1) ! ! Local: integer :: im C **** C **** CARRY OUT INTERPOLATION C **** DO 1 IM = 1,IMAXM AM(IM,1) = 1 AG(IG(IM,JM),JG(IM,JM))*WT(1,IM,JM)+ 2 AG(IG(IM,JM)+1,JG(IM,JM))*WT(2,IM,JM)+ 3 AG(IG(IM,JM)+1,JG(IM,JM)+1)*WT(3,IM,JM)+ 4 AG(IG(IM,JM),JG(IM,JM)+1)*WT(4,IM,JM) if (iprint > 0) write(6,"('grdint: i=',i3,' lat=',i3,' long=', | i3,' latg=',i3,' wght=',4e12.4,' fgeo=',e12.4,' fmag=', | e12.4)") im,jm,ig(im,jm),jg(im,jm),wt(:,im,jm), | ag(ig(im,jm),jg(im,jm)),am(im,1) 1 CONTINUE RETURN END !----------------------------------------------------------------------- SUBROUTINE ADVECL(NX,S,K) use cons_module,only: imaxp2,dphi_2div3,dphi_1div12,re_inv, | dlamda_2div3,dlamda_1div12 implicit none #include "params.h" #include "fcom.h" #include "buff.h" #include "index.h" #include "phys.h" ! ! Args: real,intent(out) :: S(1) integer,intent(in) :: nx,k ! ! Local: integer :: nuk,nxjp2k,nxjp1k,nxk,nxjm1k,nxjm2k,nvk,nvjp2k, | nvjp1k,nvjm1k,nvjm2k,i ! NUK=NJ+NU+(K-1) NXJP2K=NJP2+NX NXJP1K=NJP1+NX NXK=NJ+NX NXJM1K=NJM1+NX NXJM2K=NJM2+NX NVK=NV+K-1 NVJP2K=NJP2+NVK NVJP1K=NJP1+NVK NVJM1K=NJM1+NVK NVJM2K=NJM2+NVK DO 1 I=3,IMAXP2 S(I)=.5*RACS*(dlamda_2div3*(F(I+1,NXK)-F(I-1,NXK))*(F(I+1,NUK)+ | F(I-1,NUK))-dlamda_1div12*(F(I+2,NXK)-F(I-2,NXK))*(F(I+2,NUK)+ | F(I-2,NUK)))+.5*re_inv*(dphi_2div3*(F(I,NXJP1K)-F(I,NXJM1K))* | (F(I,NVJP1K)+F(I,NVJM1K))-dphi_1div12*(F(I,NXJP2K)-F(I,NXJM2K))* | (F(I,NVJP2K)+F(I,NVJM2K))) 1 CONTINUE RETURN END !----------------------------------------------------------------------- SUBROUTINE PRJ1(NX,RESF,NCX,RHSC,NXA,NXB) SAVE DIMENSION RHSC(NCX),RESF(NX) C C FULL WEIGHTING RESIDUAL RESTRICTION FROM RESF TO RHSC IN ONE-DIMEN C RHSC(1) = 0.0 RHSC(NCX) = 0.0 C C END POINTS C RHSC(1) = 0.0 IF (NXA.EQ.2) THEN RHSC(1) = 0.5*(RESF(1)+RESF(2)) END IF IF (NXA.EQ.0) THEN RHSC(1) = 0.25*(RESF(NX-1)+2.*RESF(1)+RESF(2)) END IF RHSC(NCX) = 0.0 IF (NXB.EQ.2) THEN RHSC(NCX) = 0.5*(RESF(NX-1)+RESF(NX)) END IF IF (NXB.EQ.0) THEN RHSC(NCX) = 0.25*(RESF(NX-1)+2.*RESF(NX)+RESF(2)) END IF C C INTERIOR C DO 1 I=3,NX-2,2 IC = I/2+1 RHSC(IC) = 0.25*(RESF(I-1)+2.*RESF(I)+RESF(I+1)) 1 CONTINUE RETURN END !----------------------------------------------------------------------- SUBROUTINE CUB1(NX,PHI,NCX,PHIC,NXA) SAVE C C ONE DIMENSIONAL CUBIC INTERPOLATION OF PHIC INTO PHI C DIMENSION PHI(NX),PHIC(NCX) COMMON/IJSF/ISTART,IFINAL,JSTART,JFINAL,ISRT,IFNL,JSRT,JFNL DO 1 I=ISTART,IFINAL,2 IC = I/2+1 PHI(I) = PHIC(IC) 1 CONTINUE C C SET EVEN X POINTS INTERPOLATING FROM ODD POINTS IN PHI C IF (NXA.NE.0) THEN C C ASYMMETRIC FORMULA NEAR NONPERIODIC BNDYS C PHI(2)=(5.*PHI(1)+15.*PHI(3)-5.*PHI(5)+PHI(7))/16. PHI(NX-1)=(5.*PHI(NX)+15.*PHI(NX-2)-5.*PHI(NX-4)+PHI(NX-6))/16. ELSE C C PERIODICITY IN X ALOWS SYMMETRIC FORMULA NEAR BNDYS C PHI(2) = (-PHI(NX-2)+9.*(PHI(1)+PHI(3))-PHI(5))/16. PHI(NX-1) = (-PHI(NX-4)+9.*(PHI(NX-2)+PHI(NX))-PHI(3))/16. END IF C C DEEP INTERIOR C DO 2 I=4,NX-3,2 PHI(I)=(-PHI(I-3)+9.*(PHI(I-1)+PHI(I+1))-PHI(I+3))/16. 2 CONTINUE RETURN END !----------------------------------------------------------------------- SUBROUTINE RES1CR(NX,NY,PHI,COF,RESF,J) SAVE C C COMPUTE RESIDUAL FOR ALL I AT CURRENT J IN RESF C COMMON/IMUD2/INTL,NXA,NXB,NYC,NYD,NXL,NYL,NSTEP,NFX,NFY,IGUESS, + MAXIT,METHOD,NWORK,LWORK,ITERO,KLEVEL,KCUR, + KCYCLE,IPRER,IPOST,IRESW COMMON/FMUD2/XA,XB,YC,YD,TOLMAX,RELMAX DIMENSION PHI(NX,NY),COF(NX,NY,10),RESF(NX) C C SET J-,J+ ALLOWING FOR PERIODIC B.C. C JM1 = J-1 IF (J.EQ.1) JM1 = NY-1 JP1 = J+1 IF (J.EQ.NY) JP1 = 2 C C SET RESIDUAL ON INTERIOR IN X C DO 1 I=2,NX-1 RESF(I) = COF(I,J,10)-( + COF(I,J,1)*PHI(I+1,J)+ + COF(I,J,2)*PHI(I+1,JP1)+ + COF(I,J,3)*PHI(I,JP1)+ + COF(I,J,4)*PHI(I-1,JP1)+ + COF(I,J,5)*PHI(I-1,J)+ + COF(I,J,6)*PHI(I-1,JM1)+ + COF(I,J,7)*PHI(I,JM1)+ + COF(I,J,8)*PHI(I+1,JM1)+ + COF(I,J,9)*PHI(I,J)) 1 CONTINUE C C SET RESIDUAL AT BOUNDARIES C RESF(1) = 0.0 RESF(NX) = 0.0 IM1 = NX-1 IP1 = 2 DO 2 I=1,NX,NX-1 IF ((I.EQ.1.AND.NXA.NE.1).OR.(I.EQ.NX.AND.NXB.NE.1)) THEN RESF(I) = COF(I,J,10)-( + COF(I,J,1)*PHI(IP1,J)+ + COF(I,J,2)*PHI(IP1,JP1)+ + COF(I,J,3)*PHI(I,JP1)+ + COF(I,J,4)*PHI(IM1,JP1)+ + COF(I,J,5)*PHI(IM1,J)+ + COF(I,J,6)*PHI(IM1,JM1)+ + COF(I,J,7)*PHI(I,JM1)+ + COF(I,J,8)*PHI(IP1,JM1)+ + COF(I,J,9)*PHI(I,J)) END IF 2 CONTINUE RETURN END