text:
!@sum  STRAT_DIAG file for special E-P flux diagnostics from strat model
!@auth B. Suozzo/J/ Lerner
!@ver  1.0

      SUBROUTINE EPFLUX (U,V,T,P)
!@sum  EPFLUX calculates finite difference EP Flux on B-grid
!@auth B. Suozzo/J. Lerner
!@ver  1.0
C**** B-grid
C**** INPUT:   U - Zonal wind (corners) (m s-1)
C****          V - Meridional wind (corners) (m s-1)
C****          W - dp/dt (downward) (P-T grid level edges) (mb m2 s-1)
C****                (L=1 is ground level)
C****          T - Potential Temperature (K)  real T = TH*(P(mb))**K
C****          P - Pressure (mb) of column from 100mb to surface
C**** OUTPUT:  FMY,FEY - North. flux of mean,eddy ang. mom.   m3 s-2
C****          FMZ,FEZ - Vert. flux of mean,eddy ang. mom.    m3 mb s-2
C****         COR,CORR - Coriolis term and transf. coriolis term  m3 s-2
C****          FER1 - North. flux of error term 1   m2 s-2
C****          ER21,ER22 - error term 2 parts 1 & 2  m s-2
C****       FMYR,FEYR - Nor. flux of transf. mean,eddy ang. mom. m3 s-2
C****     FMZR,FEZR - Vert. flx of transf. mean,eddy ang. mom. m3 mb s-2
C****          RX - <v'th'>/<dth/dp> (U-V grid level edges)  (mb m s-1)
C****                (L=1 is ground level)
C****
C**** Note: W(,,1) is really PIT, the pressure tendency, but is not used
C****
      USE MODEL_COM, only : im,jm,lm,sig,dsig,sige,psfmpt,byim
      USE GEOM, only : dxv,rapvn,rapvs,fcor,dxyv,cosv,cosp
      USE DAGCOM, only : ajl,kajl,kep,pl=>plm
      USE DYNAMICS, only : w=>conv     ! I think this is right....?

      IMPLICIT NONE
      REAL*8, INTENT(IN), DIMENSION(IM,JM,LM) :: U,V,T
      REAL*8, INTENT(IN), DIMENSION(IM,JM) :: P

C**** NOTE: AEP was a separate array but is now saved in AJL (pointer?)
c      REAL*8, DIMENSION(JM,LM,KEP) :: AEP

C**** ARRAYS CALCULATED HERE:
      REAL*8 UV(IM,JM,LM),UW(IM,JM,LM),
     *   STB(JM,LM), RXCL(LM),UCL(LM), WXXS(IM),WXXN(IM)
C**** The following quantities are added into AEP(,,1..19)
     *  ,FMY(JM,LM), FEY(JM,LM),  FMZ(JM,LM), FEZ(JM,LM)
     *  ,FMYR(JM,LM),FEYR(JM,LM), FMZR(JM,LM),FEZR(JM,LM)
     *  ,COR(JM,LM), CORR(JM,LM), FER1(JM,LM),ER21(JM,LM)
     *  ,ER22(JM,LM),VR(JM,LM),   WR(JM,LM),  RX(JM,LM)
     *  ,UI(JM,LM),  VI(JM,LM),   WI(JM,LM)
C**** End of quantities added into AEP
     *  ,TI(JM,LM),  FD(IM,JM)
C**** The following quantities are added into AEP(,,1..19)
C**** (equivalenced to XEP)
      COMMON /EPCOM1/ FMY, FEY,  FMZ, FEZ, FMYR, FEYR, FMZR, FEZR, COR,
     *     CORR, FER1, ER21, ER22, VR, WR, RX, UI, VI,  WI
      REAL*8 DUT(JM,LM), AX(JM,LM)
      REAL*8 XEP(JM,LM,KEP)
      EQUIVALENCE (XEP,FMY)

      INTEGER I,J,L,N,IM1,IP1
      REAL*8 UCB,UCT,UCM,ALPH,RXC

C     Use IDACC(4) for calling frequency

C**** Initialise for this call
      XEP = 0

C****
C**** ZONAL AVERAGE QUANTITIES
C****
C**** UI(J,L) ... ZONAL AVERAGE U-WIND (m s-1)
C**** VI(J,L) ... ZONAL AVERAGE V-WIND (m s-1)
C**** TI(J,L) ... ZONAL AVERAGE POTENTIAL TEMPERATURE (K)
C**** WI(J,L) ... ZONAL AVERAGE VERTICAL WIND (mb m2 s-1)
      CALL AVGVI (U,UI)
      CALL AVGVI (V,VI)
      CALL AVGI (W,WI)
      CALL AVGI (T,TI)
C****
C**** STB(J,L) ... DELTA THETA             (PRESSURE EDGES)
C****
      DO L=2,LM
        DO J=1,JM
          STB(J,L) = TI(J,L)-TI(J,L-1)
          IF(STB(J,L).LT.1d-2) THEN
CW         IF(STB(J,L).LT.1d-3)
CW   *     WRITE (*,*) 'STB < .001 AT J,L,STB:',J,L,STB(J,L)
            STB(J,L)=1d-2
          ENDIF
        END DO
      END DO
      STB(:,1)  = 1d-2
C****
C**** CORRELATIONS OF VARIOUS QUANTITIES
C****
C****
C**** FMY .............. (VUdx) NORTHWARD MEAN MOMENTUM TRANSPORT
C****                         (m3 s-1)
      DO L=1,LM
        FMY(1,L)=0.
        FMY(JM,L)=0.
        DO J=2,JM-1
          FMY(J,L)=.25*(DXV(J)*VI(J,L)+DXV(J+1)*VI(J+1,L))
     *         * (UI(J,L)+UI(J+1,L))
        END DO
C****
C**** FEY, UV(I,J,L) ... (V'U'dx) NORTHWARD EDDY MOMENTUM TRANSPORT
C****                           (m3 s-2)
        DO I=1,IM
          UV(I,1,L)=0.
          UV(I,JM,L)=0.
        END DO
        DO J=2,JM-1
          IM1=IM-1
          I=IM
          DO IP1=1,IM
            UV(I,J,L) = (
     *           DXV(J)*(V(IM1,J,L)+2*V(I,J,L)+V(IP1,J,L)-4*VI(J,L)) +
     *           DXV(J+1)*(V(IM1,J+1,L)+2*V(I,J+1,L)+V(IP1,J+1,L)-4*VI(J
     *           +1,L)))* (U(I,J,L)+U(I,J+1,L)-(UI(J,L)+UI(J+1,L)))
            IM1=I
            I=IP1
          END DO
        END DO
      END DO
      CALL AVGI (UV,FEY)
      DO L=1,LM
      DO J=1,JM
        FEY(J,L)=.0625d0*FEY(J,L)
      END DO
      END DO
C****
C**** FMZ .............. (WUdA) VERTICAL MEAN MOMENTUM TRANSPORT
C**** FEZ, UW(I,J,L) ... (U'W'dA) VERTICAL EDDY MOMENTUM TRANSPORT
C****                           (m3 mb s-2)
      DO L=2,LM
        WXXS(1:IM)=0.
        DO J=2,JM
          FMZ(J,L) = (WI(J-1,L)*RAPVN(J-1)+WI(J,L)*RAPVS(J))
     *         * (UI(J,L-1)+UI(J,L))
          I=IM
          DO IP1=1,IM
            WXXN(I)=(W(I,J,L)+W(IP1,J,L))-2*WI(J,L)
            UW(I,J,L) = (WXXS(I)*RAPVN(J-1)+WXXN(I)*RAPVS(J))
     *           * ((U(I,J,L-1)+U(I,J,L))-(UI(J,L-1)+UI(J,L)))
            WXXS(I)=WXXN(I)
            I=IP1
          END DO
        END DO
      END DO
      CALL AVGI (UW,FEZ)
      FMZ(:,1)=0.
      FEZ(:,1)=0.
      FMZ(1,2:LM)=0.
      FEZ(1,2:LM)=0.
      DO L=1,LM
      DO J=1,JM
        FEZ(J,L)=.5*FEZ(J,L)
      END DO
      END DO
C****
C**** COR(J,L) .......... CORIOLIS FORCE
C****                        (m3 s-2)
      DO L=1,LM
        COR( 2,L)=.5*(2.*FCOR( 1)+FCOR(   2))*VI(2,L)
        COR(JM,L)=.5*(2.*FCOR(JM)+FCOR(JM-1))*VI(JM,L)
        DO J=3,JM-1
          COR(J,L)=.5*(FCOR(J-1)+FCOR(J))*VI(J,L)
        END DO
      END DO
C****
C**** ERROR TERMS
C****
      DO L=1,LM
        FER1(:,L)=0.
        ER21(:,L)=0.
        ER22(:,L)=0.
C**** ERROR TERM 1 IS DIVERGENCE OF FER1 (m2 s-2)
        IM1=IM-1
        I=IM
        DO IP1=1,IM
          FER1(1,L) =FER1(1,L) +(U(I,2,L)**2   - U(I,2,L)**2)
          FER1(JM,L)=FER1(JM,L)+(U(IM1,JM,L)**2- U(IP1,JM,L)**2)
          DO J=2,JM-1
C? FER1 does not include COSP, yet ER1 = 1/COSV*(FER1(J)-FER1(J-1))
            FER1(J,L)= FER1(J,L) + ((U(IM1,J,L)+U(I,J+1,L))**2
     *           - (U(IP1,J,L)+U(I,J+1,L))**2)
          END DO
          IM1=I
          I=IP1
        END DO
C**** ERROR TERM 2 - THE METRIC TERM AS IN RUN 999
        IM1=IM
        DO I=1,IM
          FD(I,1)= 0.
          FD(I,JM)=0.
          DO J=2,JM-1
            FD(I,J) = (DXV(J)-DXV(J+1))*
     *           .25*(U(IM1,J,L)+U(I,J,L)+U(IM1,J+1,L)+U(I,J+1,L))
          END DO
          DO J=2,JM
            ALPH=.25*(FD(I,J-1)+FD(I,J))
            ER22(J,L)=ER22(J,L)+ALPH*(V(IM1,J,L)+V(I,J,L))
          END DO
          IM1=I
        END DO
CE*** ERROR TERM 2, PART 2  -  ER22  (m s-2)
CE    IM1=IM-1
CE    I=IM
CE    DO IP1=1,IM
CE    UXXNS=(U(IM1,2,L)+2*U(I,2,L)+U(IP1,2,L))
CE    XUXXYS=2*UXXNS*DXP(2)
CE    DO J=2,JM-1
CE    UXXNN=(U(IM1,J+1,L)+2*U(I,J+1,L)+U(IP1,J+1,L))
CE    XUXXYN=(UXXNS+UXXNN)*(DXP(J+1)-DXP(J-1))
CE    ER22(J,L) = ER22(J,L)-V(I,J,L)*(XUXXYS+XUXXYN)
CE    XUXXYS=XUXXYN
CE    UXXNS=UXXNN
CE    END DO
CE    XUXXYN=-2*UXXNS*DXP(JM-1)
CE    ER22(JM,L) = ER22(JM,L)-V(I,JM,L)*(XUXXYS+XUXXYN)
CE    IM1=I
CE    I=IP1
CE    END DO
C**** NORMALIZE ERROR TERMS AND CONVERT TO OUTPUT UNITS
      DO J=1,JM
        FER1(J,L)=1./(48*IM)*FER1(J,L)
      END DO
CE      DO J=2,JM
CE        ER21(J,L) = -1./(2*DXYV(J)*COSV(J))*
CE   *  (FMY(J-1,L)+FEY(J-1,L)+FMY(J,L)+FEY(J,L))*(COSP(J)-COSP(J-1))
CE        ER22(J,L) = 1./(32*IM*DXYV(J))*ER22(J,L)
CE      END DO
      DO J=2,JM
        ER22(J,L) = 1./(IM*DXYV(J))*ER22(J,L)
      END DO
      END DO
C****
C****  TRANSFORMED CIRCULATION
C****

C****
C**** RX(J,L) ..... ([V'TH']/[DTH/DP]) TRANSFORMATION GAUGE
C****               (U-wind grid, level edges)  (m mb s-1)
      DO L=2,LM
      DO J=2,JM
        RX(J,L) = 0.
        I=IM
        DO IP1=1,IM
          RX(J,L) = RX(J,L) + (V(I,J,L)+V(IP1,J,L)-2*VI(J,L))
     *         * (T(I,J-1,L)+T(I,J,L) - (TI(J-1,L)+TI(J,L)))
     *         +              (V(I,J,L-1)+V(IP1,J,L-1)-2*VI(J,L-1))
     *         * (T(I,J-1,L-1)+T(I,J,L-1) - (TI(J-1,L-1)+TI(J,L-1)))
          I=IP1
        END DO
        RX(J,L) = .25*BYIM * RX(J,L) *
     *       (PL(L)-PL(L-1))/(STB(J-1,L)+STB(J,L))
      END DO
      END DO
      RX(2:JM,1)=0.
C****
C**** VR(J,L) ........  TRANSFORMED WIND
C**** WR(J,L) ........
C****
      DO L=1,LM-1
      DO J=2,JM
        VR(J,L)=VI(J,L)+(RX(J,L+1)-RX(J,L))/(PSFMPT*DSIG(L))
      END DO
      END DO
      DO J=2,JM
        VR(J,LM)=VI(J,LM)-RX(J,LM)/(PSFMPT*DSIG(LM))
      END DO
      DO L=1,LM
      DO J=2,JM-1
        WR(J,L)=WI(J,L)+(RX(J+1,L)*DXV(J+1)-RX(J,L)*DXV(J))
      END DO
      END DO
      DO L=1,LM
        WR(1,L)=WI(1,L) + RX(2,L)*DXV(2)
        WR(JM,L)=WI(JM,L) - RX(JM,L)*DXV(JM)
      END DO
C****
C**** TRANSFORMED MEAN FLUXES
C****
C**** FMYR(J,L)     (m3 s-2)
C**** FMZR(J,L)     (m3 mb s-2)
C**** CORR(J,L)     (m3 s-2)
C****
      DO L=1,LM
      DO J=2,JM-1
        FMYR(J,L)=.25*(VR(J,L)*DXV(J)+VR(J+1,L)*DXV(J+1))
     *       *  (UI(J,L)+UI(J+1,L))
      END DO
      END DO
      DO L=2,LM
      DO J=2,JM
        FMZR(J,L) = (WR(J-1,L)*RAPVN(J-1)+WR(J,L)*RAPVS(J))
     *       *  (UI(J,L-1)+UI(J,L))
      END DO
      END DO
      DO L=1,LM
        CORR( 2,L)=.5*(2.*FCOR( 1)+FCOR(   2))*VR(2,L)
        CORR(JM,L)=.5*(2.*FCOR(JM)+FCOR(JM-1))*VR(JM,L)
        DO J=3,JM-1
          CORR(J,L)=.5*(FCOR(J-1)+FCOR(J))*VR(J,L)
        END DO
      END DO
C****
C**** TRANSFORMED EDDY FLUXES
C****
C**** FEYR(J,L)
C****
      DO L=1,LM
        RXCL(L)=RX(2,L)*COSV(2)
        UCL(L)=UI(2,L)*DXV(2)
      END DO
      DO J=2,JM-1
        UCB=(UI(J+1,1)*DXV(J+1)+UCL(1))
        UCM=(UI(J+1,2)*DXV(J+1)+UCL(2))
        FEYR(J,1)=0.
        FEYR(J,LM)=0.
        DO L=2,LM-1
          UCT=(UI(J+1,L+1)*DXV(J+1)+UCL(L+1))
          FEYR(J,L) = .125d0/(COSP(J)*PSFMPT*DSIG(L))
     *         *    ((RX(J+1,L)*COSV(J+1)+RXCL(L))*(UCM-UCB)
     *         +    (RX(J+1,L+1)*COSV(J+1)+RXCL(L+1))*(UCT-UCM))
          UCB=UCM
          UCM=UCT
        END DO
        DO L=1,LM
          RXCL(L)=RX(J+1,L)*COSV(J+1)
          UCL(L)=UI(J+1,L)*DXV(J+1)
        END DO
      END DO
C****
C**** FEZR
C****
      RXCL(:)=0.
      UCL(:)=0.
      DO J=2,JM-1
      DO L=2,LM
        RXC = RX(J,L)*COSV(J)
        UCB = (UI(J,L)+UI(J,L-1))*DXV(J)
        FEZR(J,L) = -.5*(FCOR(J-1)+FCOR(J))*RX(J,L)
     *       + .125d0/COSV(J)* ((RXC+RXCL(L))
     *       *     ((UI(J,L)+UI(J,L-1))*DXV(J) - (UCL(L)+UCL(L-1)))
     *       +     (RX(J+1,L)*COSV(J+1)+RXC)
     *       *     ((UI(J+1,L)+UI(J+1,L-1))*DXV(J+1) - UCB))
      END DO
      DO L=1,LM
        RXCL(L)=RX(J,L)*COSV(J)
        UCL(L)=UI(J,L)*DXV(J)
      END DO
      DO L=2,LM
        RXC=RX(JM,L)*COSV(JM)
        UCB=(UI(JM,L)+UI(JM,L-1))*DXV(JM)
        FEZR(JM,L) = -FCOR(JM)*RX(JM,L)
     *       + .125d0/COSV(JM)* (RXC+RXCL(L))
     *       *     ((UI(JM,L)+UI(JM,L-1))*DXV(JM) - (UCL(L)+UCL(L-1)))
      END DO
      END DO
C**** Add Eulerian circulation to transformed eddy components
      DO L=1,LM
      DO J=1,JM
        FEYR(J,L)=FEY(J,L)+FEYR(J,L)
        FEZR(J,L)=FEZ(J,L)+FEZR(J,L)
      END DO
      END DO
C****
C**** ACCUMULATE EP FLUXES
C****
      DO N=1,KEP-2
      DO L=1,LM
      DO J=1,JM
        AJL(J,L,KAJL-KEP+N)=AJL(J,L,KAJL-KEP+N)+XEP(J,L,N)
c        AEP(J,L,N)=AEP(J,L,N)+XEP(J,L,N)
      END DO
      END DO
      END DO
      RETURN
C****
      ENTRY EPFLXI (U)
      CALL AVGVI (U,AJL(1,1,KAJL))
c      CALL AVGVI (U,AEP(1,1,KEP))
CW          WRITE (36,'('' TAU='',F12.0)') TAU
CW          CALL WRITJL ('U - INITIAL     ',AEP(1,1,KEP),1.)
      RETURN
C****
      ENTRY EPFLXF (U)
      CALL AVGVI (U,AX)
      DO L=1,LM
      DO J=2,JM
        AJL(J,L,KAJL-1) = AX(J,L)-AJL(J,L,KAJL)
c       AEP(J,L,KEP-1) = AX(J,L)-AEP(J,L,KEP)
      END DO
      END DO
CW          CALL WRITJL ('U - FINAL       ',AX,1.)
CW          CALL WRITJL ('DU - TOTAL      ',AEP(1,1,KEP-1,1.)
CW         WRITE (36,'('' TAU='',F12.0,''  IDUM(1)='',I6)') TAU,IDUM(1)
      RETURN
      END

      SUBROUTINE EPFLXP
!@sum  EPFLXP prints out diagnostics of E-P Fluxes
!@auth B. Suozzo/J. Lerner
!@ver  1.0
C****
C**** B-grid
C**** INPUT: AEP contains the following quantities summed over time:
C****        FMY,FEY - North. flux of mean,eddy ang. mom.   m3 s-2
C****        FMZ,FEZ - Vert. flux of mean,eddy ang. mom.    m3 mb s-2
C****        COR,CORR - coriolis term and transf. coriolis term  m3 s-2
C****        FER1 - North. flux of error term 1   m2 s-2
C****        ER21,ER22 - error term 2 parts 1 & 2  m s-2
C****     FMYR,FEYR - Nor. flux of transf. mean,eddy ang. mom. m3 s-2
C****   FMZR,FEZR - Vert. flx of transf. mean,eddy ang. mom. m3 mb s-2
C**** OUTPUT:
C****        DMF,DEF - Divergence of mean,eddy ang. mom.  m3 s-2
C****        DMFR,DEFR - Div. of transf. mean,eddy ang. mom.  m3 s-2
C****        COR,CORR - same as above   m3 s-2
C****        ER1,ER2 - error terms 1 and 2   m s-2
C****   DUD,DUR - Delta U by Eulerian and transf. circulation  m s-2
C****
      USE MODEL_COM, only : im,jm,lm,sig,dsig,dtsrce=>dtsrc,psfmpt,fim
     *     ,idacc,ndaa,ls1,ptop
      USE GEOM, only : dxyv,bydxyv,cosv,cosp,dxv,dyv
      USE DAGCOM, only : ajl,kajl,kep,apj
     &     ,jl_dudfmdrg,jl_dumtndrg,jl_dushrdrg
     &     ,jl_dumcdrgm10,jl_dumcdrgp10
     &     ,jl_dumcdrgm40,jl_dumcdrgp40
     &     ,jl_dumcdrgm20,jl_dumcdrgp20
     &     ,jl_dudtsdif,jl_damdc,jl_dammc
      IMPLICIT NONE
C**** NOTE: AEP was a separate array but is now saved in AJL (pointer?)
c      REAL*8, DIMENSION(JM,LM,KEP) :: AEP

c      COMMON /PROGCB/ U,V,T,SX,SY,SZ,P,Q   !not used?
C**** diagnostic information for print out
C**** this should be in an init_ep routine or something
      integer, parameter :: njl_out=7
      character(len=50) :: lname(njl_out) = (/
     *     'DU/DT BY EULER CIRC. + CONVEC + DRAG+DIF+ER2',
     *     'DU/DT BY MEAN ADVECTION                     ',
     *     'DU/DT BY EDDY CONVERGENCE                   ',
     *     'DU/DT BY TRANSFORMED ADVECTION              ',
     *     'DU/DT BY ELIASSEN-PALM DIVERGENCE           ',
     *     'DU/DT BY F.D. ERROR TERM 1                  ',
     *     'DU/DT BY F.D. ERROR TERM 2                  '/)
      character(len=30) :: sname(njl_out) = (/
     *     'tot_dudt     ','dudt_meanadv ','dudt_eddycnv '
     *     ,'dudt_trnsadv ','dudt_epflxdiv','dudt_fderr1  '
     *     ,'dudt_fderr2  '/)
      character(len=50) :: units(njl_out) = 'm/s^2'
      integer, dimension(njl_out) :: pow = (/ -6,-6,-6,-6,-6,-6,-6 /)

C**** ARRAYS CALCULATED HERE:
      REAL*8 DXCOSV(JM), ONES(JM+LM),PMO(LM),DP(LM)
     *  ,FMY(JM,LM), FEY(JM,LM),  FMZ(JM,LM), FEZ(JM,LM)
     *  ,FMYR(JM,LM),FEYR(JM,LM), FMZR(JM,LM),FEZR(JM,LM)
     *  ,COR(JM,LM), CORR(JM,LM), FER1(JM,LM),ER21(JM,LM)
     *  ,ER22(JM,LM),VR(JM,LM),   WR(JM,LM),  RX(JM,LM)
     *  ,UI(JM,LM),  VI(JM,LM),   WI(JM,LM),  DUT(JM,LM)
     *  ,DUD(JM,LM), DUDS(JM,LM), DUR(JM,LM)
     *  ,DMF(JM,LM),DEF(JM,LM),DMFR(JM,LM),DEFR(JM,LM)
CW   *  ,DMY(JM,LM),DMZ(JM,LM),DEY(JM,LM),DEZ(JM,LM)
CW     *  ,DMYR(JM,LM),DMZR(JM,LM),DEYR(JM,LM),DEZR(JM,LM)
     *  ,ER1(JM,LM),ER2(JM,LM), BYDPJL(JM,LM)
C**** common needed to pass from XEP to individual arrays (better way??)
C**** Not clear how many of these are actually used.
      COMMON /EPCOM/ FMY, FEY, FMZ, FEZ, FMYR, FEYR, FMZR, FEZR, COR,
     *     CORR, FER1, ER21, ER22, VR, WR, RX, UI, VI, WI, DUT, DUD
      REAL*8, DIMENSION(JM,LM,KEP) :: XEP
      EQUIVALENCE (XEP,FMY)
      REAL*8 DTAEP,BYDT,SCALEP,SCALE1,BYIAEP
      INTEGER I,J,L,N,JL

C**** Initialize constants
      DTAEP = DTsrce*NDAA   ! change of definition of NDAA
      BYIAEP=1./(IDACC(4)+1.D-20)
      BYDT=1./(DTSRCE*IDACC(1)+1.D-20)
      DO J=2,JM
        DXCOSV(J) = DXV(J)*COSV(J)
      END DO
      DO JL=1,JM+LM
        ONES(JL)=1.
      END DO
      DO L=1,LM
        DP(L)  = PSFMPT*DSIG(L)
        PMO(L) = PSFMPT*SIG(L)+PTOP
      END DO
      DO L=1,LS1-1
      DO J=2,JM
        BYDPJL(J,L)=(.25*FIM*IDACC(4))/(APJ(J,2)+1.D-20)
      END DO
      END DO
      DO L=LS1,LM
      DO J=2,JM
        BYDPJL(J,L)=1./DP(L)
      END DO
      END DO
C**** Normalize AEP
C     CALL EPFLXF (U)
      DO N=1,KEP-2
      DO L=1,LM
      DO J=1,JM
        XEP(J,L,N)=AJL(J,L,KAJL-KEP+N)*BYIAEP
c        XEP(J,L,N)=AEP(J,L,N)*BYIAEP
      END DO
      END DO
      END DO
      DO N=KEP-1,KEP
      DO L=1,LM
      DO J=1,JM
        XEP(J,L,N)=AJL(J,L,KAJL-KEP+N)
c        XEP(J,L,N)=AEP(J,L,N)
      END DO
      END DO
      END DO
C****
C**** DMF,DEF by horizontal convergence
C****
      DO L=1,LM
      DO J=2,JM
        DMF(J,L)  = (FMY(J-1,L)*COSP(J-1)-FMY(J,L)*COSP(J))/COSV(J)
        DEF(J,L)  = (FEY(J-1,L)*COSP(J-1)-FEY(J,L)*COSP(J))/COSV(J)
C       DMY(J,L)  = (FMY(J-1,L)*COSP(J-1)-FMY(J,L)*COSP(J))/COSV(J)
C       DEY(J,L)  = (FEY(J-1,L)*COSP(J-1)-FEY(J,L)*COSP(J))/COSV(J)
C       DMYR(J,L) = (FMYR(J-1,L)*COSP(J-1)-FMYR(J,L)*COSP(J))/COSV(J)
C       DEYR(J,L) = (FEYR(J-1,L)*COSP(J-1)-FEYR(J,L)*COSP(J))/COSV(J)
        DMFR(J,L) = (FMYR(J-1,L)*COSP(J-1)-FMYR(J,L)*COSP(J))/COSV(J)
        DEFR(J,L) = (FEYR(J-1,L)*COSP(J-1)-FEYR(J,L)*COSP(J))/COSV(J)
      END DO
      END DO
C****
C**** Add DMF,DEF by vertical convergence
C****
      DO J=2,JM
        DMF(J,LM)  = DMF(J,LM)  - FMZ(J,LM)/BYDPJL(J,LM)
        DEF(J,LM)  = DEF(J,LM)  - FEZ(J,LM)/BYDPJL(J,LM)
CW      DMZ(J,LM)  = - FMZ(J,LM)/BYDPJL(J,LM)
CW      DEZ(J,LM)  = - FEZ(J,LM)/BYDPJL(J,LM)
CW      DMZR(J,LM) = - FMZR(J,LM)/BYDPJL(J,LM)
CW      DEZR(J,LM) = - FEZR(J,LM)/BYDPJL(J,LM)
        DMFR(J,LM) = DMFR(J,LM) - FMZR(J,LM)/BYDPJL(J,LM)
        DEFR(J,LM) = DEFR(J,LM) - FEZR(J,LM)/BYDPJL(J,LM)
        DO L=1,LM-1
          DMF(J,L)  = DMF(J,L)  + (FMZ(J,L+1)-FMZ(J,L))*BYDPJL(J,L)
          DEF(J,L)  = DEF(J,L)  + (FEZ(J,L+1)-FEZ(J,L))*BYDPJL(J,L)
CW        DMZ(J,L)  = (FMZ(J,L+1)-FMZ(J,L))*BYDPJL(J,L)
CW        DEZ(J,L)  = (FEZ(J,L+1)-FEZ(J,L))*BYDPJL(J,L)
CW        DMZR(J,L) = (FMZR(J,L+1)-FMZR(J,L))*BYDPJL(J,L)
CW        DEZR(J,L) = (FEZR(J,L+1)-FEZR(J,L))*BYDPJL(J,L)
          DMFR(J,L) = DMFR(J,L) + (FMZR(J,L+1)-FMZR(J,L))*BYDPJL(J,L)
          DEFR(J,L) = DEFR(J,L) + (FEZR(J,L+1)-FEZR(J,L))*BYDPJL(J,L)
        END DO
      END DO
C****
C**** ADD COR,CORR
C****
      DO L=1,LM
        DO J=2,JM
          DMF(J,L)  = DMF(J,L)  + COR(J,L)
C         DMY(J,L)  = DMY(J,L)  + COR(J,L)
C         DMYR(J,L) = DMYR(J,L) + CORR(J,L)
          DMFR(J,L) = DMFR(J,L) + CORR(J,L)
        END DO
      END DO
C****
C**** ER1 by horizontal convergence (m s-2)
C****
      DO L=1,LM
      DO J=2,JM
        ER1(J,L) = 1./(DYV(J)*COSV(J))*(FER1(J-1,L)-FER1(J,L))
        ER2(J,L) = ER21(J,L)+ER22(J,L)
      END DO
      END DO
C****
C**** DUD,DUR total change in Eulerian, Transformed wind  (m s-2)
C****
      DO L=1,LM
      DO J=2,JM
        DUR(J,L)  = BYDXYV(J)*(DMFR(J,L) +DEFR(J,L))
        DUD(J,L)  = BYDXYV(J)*(DMF(J,L)+DEF(J,L))
      END DO
      END DO
C****
C**** Let Transformed == Transformed - Eulerian
C****
CW      DO L=1,LM
CW      DO J=1,JM
CW        DUR(J,L) = DUR(J,L) - DUD(J,L)
CW        DMFR(J,L) = DMFR(J,L) - DMF(J,L)
CW        DEFR(J,L) = DEFR(J,L) - DEF(J,L)
CW      END DO
CW      END DO
C****
C**** Print maps of EP fluxes
C**** note: JLMAP (lname,sname,units,power,Pres,Array,ScalP,ScalJ,ScalL)
C****          prints maps of  (Array * SCALEP * ScaleJ * ScaleL)
C****
      DO L=1,LM
      DO J=2,JM
        DUDS(J,L)=((AJL(J,L,JL_DUDFMDRG)+AJL(J,L,JL_DUMTNDRG))+
     &             (AJL(J,L,JL_DUSHRDRG)+AJL(J,L,JL_DUMCDRGM10))+
     *       ((AJL(J,L,JL_DUMCDRGP10)+AJL(J,L,JL_DUMCDRGM40))+
     &        (AJL(J,L,JL_DUMCDRGP40)+AJL(J,L,JL_DUMCDRGM20)))+
     *       (AJL(J,L,JL_DUMCDRGP20)+AJL(J,L,JL_DUDTSDIF)))
      END DO
      END DO
      SCALE1=1./(FIM*DTSRCE*IDACC(1)+1.D-20)
      DO L=1,LM
      DO J=2,JM
        DUDS(J,L) = DUDS(J,L)+
     &        (AJL(J,L,JL_DAMDC)+AJL(J,L,JL_DAMMC))*BYDPJL(J,L)
      END DO
      END DO
      DO L=1,LM
      DO J=2,JM
        DUDS(J,L) = (DUD(J,L)-ER2(J,L)) + DUDS(J,L)*SCALE1
      END DO
      END DO
      SCALEP=1
CW      CALL WRITJL ('DUDT: EUL+SOURCE',DUDS,SCALEP)
CW      CALL WRITJL ('DUDT: ENTIRE GCM',DUT,SCALEP) ! AJK-47 DIAGJK
      CALL JLMAP (LNAME(1),SNAME(1),UNITS(1),POW(1),PMO,DUDS,SCALEP,ONES
     *     ,ONES,LM,2,2)
C****
      CALL JLMAP (LNAME(2),SNAME(2),UNITS(2),POW(2),PMO,DMF,SCALEP
     *     ,BYDXYV,ONES,LM,2,2)
      CALL JLMAP (LNAME(3),SNAME(3),UNITS(3),POW(3),PMO,DEF,SCALEP
     *     ,BYDXYV,ONES,LM,2,2)
      CALL JLMAP (LNAME(4),SNAME(4),UNITS(4),POW(4),PMO,DMFR,SCALEP
     *     ,BYDXYV,ONES,LM,2,2)
      CALL JLMAP (LNAME(5),SNAME(5),UNITS(5),POW(5),PMO,DEFR,SCALEP
     *     ,BYDXYV,ONES,LM,2,2)
      CALL JLMAP (LNAME(6),SNAME(6),UNITS(6),POW(6),PMO,ER1,SCALEP,ONES
     *     ,ONES,LM,2,2)
      CALL JLMAP (LNAME(7),SNAME(7),UNITS(7),POW(7),PMO,ER2,SCALEP,ONES
     *     ,ONES,LM,2,2)
C****
CW      DO L=1,LM
CW      DO J=2,JM
CW      DMF(J,L)=DMF(J,L)*BYDXYV(J)
CW      DEF(J,L)=DEF(J,L)*BYDXYV(J)
CW      DMFR(J,L)=DMFR(J,L)*BYDXYV(J)
CW      DEFR(J,L)=DEFR(J,L)*BYDXYV(J)
CW      END DO
CW      END DO
CW      WRITE (36,'(''DU/DT = 10**-6 M S-2'')')
CW      WRITE (36,'(''TR: TRANSFORMED - EULERIAN'')')
CW      SCALEP = 1.E6
CW      CALL WRITJL ('DUDT: MEAN EULER',DMF,SCALEP)
CW      CALL WRITJL ('DUDT: EDDY EULER',DEF,SCALEP)
CW      CALL WRITJL ('DUDT: EULER CIRC',DUD,SCALEP)
CW      CALL WRITJL ('DUDT: MEAN TRANS',DMFR,SCALEP)
CW      CALL WRITJL ('DUDT: EDDY TRANS',DEFR,SCALEP)
CW      CALL WRITJL ('DUDT: TRANS-EULE',DUR,SCALEP)
      RETURN
      END

      SUBROUTINE AVGI (X,XI)
!@sum  AVGI average a 3-dimensional array in the x-direction
!@auth B. Suozzo
!@ver  1.0
      USE MODEL_COM, only : im,jm,lm,BYIM
      USE GEOM, only : imaxj
      IMPLICIT NONE

!@var X input 3-D array
      REAL*8, INTENT(IN), DIMENSION(IM,JM,LM) :: X
!@var XI output zonally averaged 2-D array
      REAL*8, INTENT(OUT), DIMENSION(JM,LM) :: XI
      INTEGER I,J,L
      REAL*8 XXI

      DO L=1,LM
        DO J=1,JM
          XXI=0.
          DO I=1,IMAXJ(J)
            XXI = XXI + X(I,J,L)
          END DO
          XI(J,L) = XXI/IMAXJ(J)
        END DO
      END DO
      RETURN
C****
      ENTRY AVGVI (X,XI)
!@sum  AVGVI average a 3-dimensional array in the x-direction (no pole)
!@auth B. Suozzo
!@ver  1.0
C****
      DO L=1,LM
        DO J=2,JM
          XXI=0.
          DO I=1,IM
            XXI = XXI + X(I,J,L)
          END DO
          XI(J,L) = XXI*BYIM
        END DO
      END DO
C****
      RETURN
      END


text:
!*******************************************************************************
!>  @file siesta_currents.f90
!>  @brief Contains routines to compute the current density from Curl(B).
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!>  Takes the Curl(B) to return the contravariant current density. Contravariant
!>  components of the magnetic field must be converted to covariant components
!>  first.
!*******************************************************************************
      MODULE siesta_currents
      USE v3_utilities, ONLY: assert
      USE stel_kinds
      USE stel_constants
      USE descriptor_mod, ONLY: SIESTA_COMM
      USE shared_data, ONLY: l_update_state, lasym, ste,                       &
                             jsju_ratio_s, jsju_ratio_a
      USE fourier
      USE nscalingtools, ONLY: nranks, startglobrow, endglobrow, MPI_ERR
      USE metrics, ONLY: tolowerf, tolowerh
      USE quantities

      IMPLICIT NONE

      CONTAINS

!*******************************************************************************
!  UTILITY SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Compute current density from Curl(B).
!>
!>  Current density K is defined as
!>
!>     K = Curl(B)/mu0                                                       (1)
!>
!>  In curvilinear coordinates this becomes
!>
!>     K^s = 1/sqrt(g)(dB_v/du - dB_u/dv)                                   (2a)
!>
!>     K^u = 1/sqrt(g)(dB_s/dv - dB_v/ds)                                   (2b)
!>
!>     K^v = 1/sqrt(g)(dB_u/ds - dB_s/du)                                   (2c)
!>
!>  where sqrt(g) is the coordinate system Jacobian. Since there are derivatives
!>  with respect to s, these routines take the half grid values as input. The
!>  contravariant current computed will be on the full grid. These routines also
!>  output currents in real space. Note the computed Kmn are [sqrt(g)*K]mn with
!>  the embedded jacobian factor.
!>
!>  @param[in]    bsupsijh B^s on the half grid.
!>  @param[in]    bsupuijh B^u on the half grid.
!>  @param[in]    bsupvijh B^v on the half grid.
!>  @param[inout] ksupsijf K^s on the full grid.
!>  @param[inout] ksupuijf K^u on the full grid.
!>  @param[inout] ksupvijf K^v on the full grid.
!>  @param[in]    lmagen   Logical controlling calc of magnetic energy wb.
!>  @param[in]    lcurr    UNKNOWN
!-------------------------------------------------------------------------------
      SUBROUTINE cv_currents(bsupsijh, bsupuijh, bsupvijh,                     &
                             ksupsijf, ksupuijf, ksupvijf,                     &
                             lmagen, lcurr)
      USE island_params, ONLY: nu_i, rmajor_i, fourier_context
      USE Hessian, ONLY: mupar_norm
      USE descriptor_mod, ONLY: iam
      USE quantities, ONLY: bsq
      USE timer_mod
      USE mpi_inc
      USE descriptor_mod, ONLY: SIESTA_COMM
      USE shared_data, ONLY: siesta_curtor, unit_out

      IMPLICIT NONE

!  Declare Arguments
      REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in)    :: bsupsijh
      REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in)    :: bsupuijh
      REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in)    :: bsupvijh
      REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupsijf
      REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupuijf
      REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupvijf
      LOGICAL, INTENT(in)                                     :: lmagen
      LOGICAL, INTENT(in)                                     :: lcurr

!  Local Variables
      INTEGER                                                 :: istat
      REAL (dp)                                               :: beta
      REAL (dp)                                               :: ton
      REAL (dp)                                               :: toff
      REAL (dp), DIMENSION(6)                                 :: temp(6)
      REAL (dp)                                               :: tmp
      INTEGER                                                 :: nsmin
      INTEGER                                                 :: nsmax
      INTEGER                                                 :: nmin
      INTEGER                                                 :: nmax
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bsubsijh
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bsubuijh
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bsubvijh
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bs_filter
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bu_filter
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bv_filter
      REAL (dp)                                               :: edge_cur

!  Start of executable code
      CALL second0(ton)

      nmin = MAX(1, startglobrow)
      nmax = MIN(ns, endglobrow)
      nsmin = LBOUND(bsupsijh,3)
      nsmax = UBOUND(bsupsijh,3)

!  Calculate covariant BsubXh on half mesh
      ALLOCATE(bsubsijh(ntheta,nzeta,nsmin:nsmax),                             &
               bsubuijh(ntheta,nzeta,nsmin:nsmax),                             &
               bsubvijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
      CALL assert_eq(0, istat, 'Allocation1 failed in cv_currents')

      CALL tolowerh(bsupsijh, bsupuijh, bsupvijh,                              &
                    bsubsijh, bsubuijh, bsubvijh, nsmin, nsmax)

!  Compute parallel damping scaling cofficient <bsupv^2/B.B> ~ k^2
      IF (.not.ALLOCATED(mupar_norm)) THEN
         ALLOCATE(mupar_norm(nsmin:nsmax), stat=istat)
         CALL assert_eq(0, istat,'mupar_norm allocation failed!')
         beta = MAX(1.E-5_dp, wp_i/wb_i)
         mupar_norm = signjac*beta/(rmajor_i*rmajor_i)  !/ vp_f(nsmin:nsmax)     !add 1/jac correction factor in AddDamping
         IF (iam .eq. 0 .and. lverbose) THEN
            WRITE (*,1000) rmajor_i
         END IF
      END IF

      IF (lmagen) THEN
         bsq(:,:,nmin:nmax) = bsupsijh(:,:,nmin:nmax)*bsubsijh(:,:,nmin:nmax)  &
                            + bsupuijh(:,:,nmin:nmax)*bsubuijh(:,:,nmin:nmax)  &
                            + bsupvijh(:,:,nmin:nmax)*bsubvijh(:,:,nmin:nmax)
         IF (nmin .EQ. 1) THEN
            bsq(:,:,1) = 0
         END IF

         CALL assert(nmin                   .eq. 1 .or.                        &
                     ALL(bsq(:,:,nmin:nmax) .gt. zero),                        &
                     'BSQ <= 0 in cv_currents!')

         wb = SUM(bsq(:,:,nmin:nmax)*jacobh(:,:,nmin:nmax)*wint(:,:,nmin:nmax))

         IF (PARSOLVER) THEN
#if defined(MPI_OPT)
!  FIXME: All reduce is not deterministic. This causes a divergent run sequence.
            CALL MPI_ALLREDUCE(MPI_IN_PLACE, wb, 1, MPI_REAL8, MPI_SUM,        &
                               SIESTA_COMM, MPI_ERR)
#endif
         END IF
         wb = signjac*(twopi*pi*hs_i)*wb
      END IF

!  Calculate BsubXmnh and KsupXmnf harmonics. Deallocated in calling routine.
      IF (l_update_state) THEN
         ALLOCATE(bs_filter(ntheta,nzeta,nsmin:nsmax),                         &
                  bu_filter(ntheta,nzeta,nsmin:nsmax),                         &
                  bv_filter(ntheta,nzeta,nsmin:nsmax), stat=istat)
         CALL assert_eq(0, istat, 'Allocation2 failed in CURLB')
      END IF

      edge_cur = 0.0
      CALL curlb(ksupsmnsf, ksupumncf, ksupvmncf,                              &
                 bsubsijh, bsubuijh, bsubvijh,                                 &
                 bs_filter, bu_filter, bv_filter, edge_cur,                    &
                 jsju_ratio_s, f_sin, lcurr, nsmax, nsmin, siesta_curtor)
      IF (lasym) THEN
         CALL curlb(ksupsmncf, ksupumnsf, ksupvmnsf,                           &
                    bsubsijh, bsubuijh, bsubvijh,                              &
                    bs_filter, bu_filter, bv_filter, edge_cur,                 &
                    jsju_ratio_a, f_cos, lcurr, nsmax, nsmin,                  &
                    siesta_curtor)
      END IF

      IF (nsmax .eq. ns) THEN
         siesta_curtor = -fourier_context%orthonorm(m0,n0)*siesta_curtor       &
                       / b_factor
      END IF

      IF (lcurr) THEN
         nsmax = MIN(ns, endglobrow + 1)
      ELSE
         nsmax = MIN(ns, endglobrow)
      END IF

!  Calculate full mesh, contravariant real-space current components
!  KsupXF = sqrt(g)*JsupXF.
      CALL getcv(ksupsmnsf, ksupumncf, ksupvmncf,                              &
                 ksupsijf, ksupuijf, ksupvijf, f_sin, nsmin, nsmax)
      IF (lasym) THEN
         CALL getcv(ksupsmncf, ksupumnsf, ksupvmnsf,                           &
                    ksupsijf, ksupuijf, ksupvijf, f_cos, nsmin, nsmax)
      END IF

!  Covariant K components (with jacob factor still there) in real space needed
!  for current diffusion
      IF (lcurr) THEN
         CALL tolowerf(ksupsijf, ksupuijf, ksupvijf,                           &
                       ksubsijf, ksubuijf, ksubvijf, nsmin, nsmax)
         IF (nsmin .eq. 1) THEN
            ksubuijf(:,:,1) = 0
         END IF
      END IF

!  Diagnostic output (only compute when state is updated)
      DIAGNO: IF (l_update_state) THEN

!  Compute spectral truncation error measure (Use ksupsijf for temporary
!  filtered bsubXijh values)
         IF (nsmax .eq. ns) THEN
            tmp = SUM(wint(:,:,ns)*(bsubvijh(:,:,ns)*bsubvijh(:,:,ns) +        &
                                    bsubuijh(:,:,ns)*bsubuijh(:,:,ns)))
            IF (tmp .ne. zero) THEN
               tmp = SQRT(edge_cur/tmp)
            END IF
         END IF
#if defined(MPI_OPT)
         IF (PARSOLVER) THEN
            CALL MPI_BCAST(tmp, 1, MPI_REAL8, nranks - 1, SIESTA_COMM, MPI_ERR)
         END IF
#endif
         IF (iam .EQ. 0) THEN
            IF (lverbose) THEN
               WRITE (*, 1001) tmp
            END IF
            WRITE (unit_out, 1001) tmp
            CALL FLUSH(unit_out)
         END IF

         temp(1) = SUM(wint(:,:,nmin:nmax) *                                   &
                       (bs_filter - bsubsijh(:,:,nmin:nmax))**2)
         temp(2) = SUM(wint(:,:,nmin:nmax) *                                   &
                       (bs_filter + bsubsijh(:,:,nmin:nmax))**2)
         temp(3) = SUM(wint(:,:,nmin:nmax) *                                   &
                       (bu_filter - bsubuijh(:,:,nmin:nmax))**2)
         temp(4) = SUM(wint(:,:,nmin:nmax) *                                   &
                       (bu_filter + bsubuijh(:,:,nmin:nmax))**2)
         temp(5) = SUM(wint(:,:,nmin:nmax) *                                   &
                       (bv_filter - bsubvijh(:,:,nmin:nmax))**2)
         temp(6) = SUM(wint(:,:,nmin:nmax) *                                   &
                       (bv_filter + bsubvijh(:,:,nmin:nmax))**2)
#if defined(MPI_OPT)
         IF (PARSOLVER) THEN
            CALL MPI_ALLREDUCE(MPI_IN_PLACE, temp, 6, MPI_REAL8, MPI_SUM,      &
                               SIESTA_COMM, MPI_ERR)
         END IF
#endif
         IF (temp(2) .ne. zero) THEN
            ste(2) = SQRT(temp(1)/temp(2))
         END IF
         IF (temp(4) .ne. zero) THEN
            ste(3) = SQRT(temp(3)/temp(4))
         END IF
         IF (temp(6) .ne. zero) THEN
            ste(4) = SQRT(temp(5)/temp(6))
         END IF
         DEALLOCATE (bs_filter, bu_filter, bv_filter)
      END IF DIAGNO

      DEALLOCATE (bsubsijh, bsubuijh, bsubvijh)

      CALL second0(toff)
      cv_current_time = cv_current_time + (toff - ton)
      time_current = time_current + (toff - ton)

1000  FORMAT(' Rmajor_i = ',1p,e12.3)
1001  FORMAT(' RMS EDGE CURRENT ||KSUPS(NS)|| = ',1p,e10.3)

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Compute fourier components of contravariant currents on full radial
!>  grid.
!>
!>  @param[inout] ksupsmnf   Contravarant s component of the current with
!>                           jacobian.
!>  @param[inout] ksupumnf   Contravarant u component of the current with
!>                           jacobian.
!>  @param[inout] ksupvmnf   Contravarant v component of the current with
!>                           jacobian.
!>  @param[in]    bsubsijh   Covariant s component of the magnetic field.
!>  @param[in]    bsubuijh   Covariant u component of the magnetic field.
!>  @param[inout] bsubvijh   Covariant v component of the magnetic field.
!>  @param[inout] bs_filter
!>  @param[inout] bu_filter
!>  @param[inout] bv_filter
!>  @param[inout] edge_cur   Edge current.
!>  @param[out]   jsju_ratio Ratio of <K^s + r12*K^u>/<K^s - r12*K^u>.
!>  @param[in]    iparity    Parity flag for the fouier components.
!>  @param[in]    lcurr      UNKNOWN
!>  @param[in]    nsmax      Maximum radial index.
!>  @param[in]    nsmin      Minimum radial index.
!>  @param[inout] curtor     Current enclosed at the boundary.
!-------------------------------------------------------------------------------
      SUBROUTINE curlb(ksupsmnf, ksupumnf, ksupvmnf,                           &
                       bsubsijh, bsubuijh, bsubvijh,                           &
                       bs_filter, bu_filter, bv_filter, edge_cur,              &
                       jsju_ratio, iparity, lcurr, nsmax, nsmin, curtor)
      USE utilities, ONLY: curl_htof
      USE island_params, ONLY: fourier_context

      IMPLICIT NONE

!  Declare Arguments
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: ksupsmnf
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: ksupumnf
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: ksupvmnf
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(in)    :: bsubsijh
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(in)    :: bsubuijh
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: bsubvijh
      REAL (dp), DIMENSION(:,:,:), INTENT(inout)              :: bs_filter
      REAL (dp), DIMENSION(:,:,:), INTENT(inout)              :: bu_filter
      REAL (dp), DIMENSION(:,:,:), INTENT(inout)              :: bv_filter
      REAL (dp), INTENT(inout)                                :: edge_cur
      REAL (dp), INTENT(out)                                  :: jsju_ratio
      INTEGER, INTENT(in)                                     :: iparity
      LOGICAL, INTENT(in)                                     :: lcurr
      INTEGER, INTENT(in)                                     :: nsmax
      INTEGER, INTENT(in)                                     :: nsmin
      REAL (dp), INTENT(inout)                                :: curtor

!  local variables
      INTEGER                                                 :: fours
      INTEGER                                                 :: fouruv
      INTEGER                                                 :: fcomb
      INTEGER                                                 :: nmin
      INTEGER                                                 :: nmax
      INTEGER                                                 :: istat
      REAL (dp)                                               :: r12
      REAL (dp), DIMENSION(-ntor:ntor)                        :: temp
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bsubsmnh
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bsubumnh
      REAL (dp), ALLOCATABLE, DIMENSION(:,:,:)                :: bsubvmnh

!  Start of executable code
      nmin = MAX(1, startglobrow)
      nmax = MIN(ns, endglobrow)

      IF (iparity .EQ. f_sin) THEN
         fcomb = f_none
         fours = f_sin
         fouruv = f_cos
         r12 = 0.5*hs_i
      ELSE
         fcomb = f_sum
         fours = f_cos
         fouruv = f_sin
         r12 = -0.5*hs_i
      END IF

      ALLOCATE(bsubsmnh(0:mpol,-ntor:ntor,nsmin:nsmax),                        &
               bsubumnh(0:mpol,-ntor:ntor,nsmin:nsmax),                        &
               bsubvmnh(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
      CALL assert_eq(0, istat, 'Allocation1 failed in CURLB')

!  Compute bsubxmnh
      CALL fourier_context%tomnsp(bsubsijh, bsubsmnh, fours)
      CALL fourier_context%tomnsp(bsubuijh, bsubumnh, fouruv)
      CALL fourier_context%tomnsp(bsubvijh, bsubvmnh, fouruv)

      IF (lcurr) THEN
         nmax = MIN(endglobrow + 1, ns)
      ELSE
         nmax = MIN(endglobrow, ns)
      END IF

      CALL CURL_HtoF(bsubsmnh, bsubumnh, bsubvmnh,                             &
                     ksupsmnf, ksupumnf, ksupvmnf,                             &
                     iparity, nsmin, nsmax, nmax, curtor)

!  Diagnostics
      DIAGNO: IF (l_update_state) THEN
         IF (nmax .eq. ns) THEN
            edge_cur = edge_cur + SUM(ksupsmnf(:,:,ns)*ksupsmnf(:,:,ns))
         END IF

         CALL fourier_context%toijsp(bsubsmnh(:,:,nmin:nmax), bs_filter,       &
                                     fcomb, fours)
         CALL fourier_context%toijsp(bsubumnh(:,:,nmin:nmax), bu_filter,       &
                                     fcomb, fouruv)
         CALL fourier_context%toijsp(bsubvmnh(:,:,nmin:nmax), bv_filter,       &
                                     fcomb, fouruv)
	     IF (nmin .eq. 1) THEN
            bv_filter(:,:,1) = bv_filter(:,:,2)
            bsubvijh(:,:,1) = bsubvijh(:,:,2)
         END IF

!  K^s-r12*K^u at origin BEFORE bc applied
         IF (nsmin .eq. 1) THEN
            temp = ksupsmnf(m1,:,1) + r12*ksupumnf(m1,:,1)
            jsju_ratio = SUM(temp*temp)
            IF (jsju_ratio .gt. zero) THEN
               temp = ksupsmnf(m1,:,1) - r12*ksupumnf(m1,:,1)
               jsju_ratio = SUM(temp*temp)/jsju_ratio
               jsju_ratio = SQRT(jsju_ratio)
            END IF
         END IF
      END IF DIAGNO

      DEALLOCATE (bsubsmnh, bsubumnh, bsubvmnh, stat=istat)
      CALL assert_eq(0, istat, 'Deallocation failed in CURLB')

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Transform full grid contravariant currents to real space.
!>
!>  @param[in]    ksupsmnf Contravariant fourier current in the s direction.
!>  @param[in]    ksupumnf Contravariant fourier current in the u direction.
!>  @param[in]    ksupvmnf Contravariant fourier current in the v direction.
!>  @param[inout] ksupsmnf Contravariant real space current in the s direction.
!>  @param[inout] ksupumnf Contravariant real space current in the u direction.
!>  @param[inout] ksupvmnf Contravariant real space current in the v direction.
!>  @param[in]    parity   Partity flag fourier components.
!>  @param[in]    nsmin    Minimum radial index.
!>  @param[in]    nsmax    Maximum radial index.
!-------------------------------------------------------------------------------
      SUBROUTINE getcv(ksupsmnf, ksupumnf, ksupvmnf,                           &
                       ksupsijf, ksupuijf, ksupvijf, parity, nsmin, nsmax)
      USE island_params, ONLY: fourier_context

      IMPLICIT NONE

!  Declare Arguments
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in)    :: ksupsmnf
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in)    :: ksupumnf
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in)    :: ksupvmnf
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupsijf
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupuijf
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupvijf
      INTEGER, INTENT(in)                                    :: parity
      INTEGER, INTENT(in)                                    :: nsmin
      INTEGER, INTENT(in)                                    :: nsmax

!  local variables
      INTEGER                                                :: fours
      INTEGER                                                :: fouruv
      INTEGER                                                :: fcomb

!  Start of executable code
      IF (parity .eq. f_sin) THEN
         fcomb = f_none
         fours = f_sin
         fouruv = f_cos
      ELSE
         fcomb = f_sum
         fours = f_cos
         fouruv = f_sin
      END IF

      CALL fourier_context%toijsp(ksupsmnf(:,:,nsmin:nsmax), ksupsijf, fcomb,  &
                                  fours)
      CALL fourier_context%toijsp(ksupumnf(:,:,nsmin:nsmax), ksupuijf, fcomb,  &
                                  fouruv)
      CALL fourier_context%toijsp(ksupvmnf(:,:,nsmin:nsmax), ksupvijf, fcomb,  &
                                  fouruv)

      END SUBROUTINE

      END MODULE

max_stars_repo_name: ORNL-Fusion/SIESTA

max_issues_repo_name: ORNL-Fusion/SIESTA

max_forks_repo_name: ORNL-Fusion/SIESTA


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module atmos_carbon_aerosol_mod
! <CONTACT EMAIL="wfc@gfdl.noaa.gov">
!   William Cooke
! </CONTACT>

! <REVIEWER EMAIL="lwh@gfdl.noaa.gov">
!   Larry Horowitz
! </REVIEWER>

! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>

! <OVERVIEW>
!     This code allows the implementation of black and organic carbon
!     tracers in the FMS framework.
! </OVERVIEW>

! <DESCRIPTION>
!   This module presents the method of Cooke et al. (1999, 2002)

!   In its present implementation the black and organic carbon tracers
!   are from the combustion of fossil fuel.

!   While the code here should provide insights into the carbonaceous
!   aerosol cycle it is provided here more as an example of how to implement
!   a tracer module in the FMS infrastructure. The parameters of the model
!   should be checked and set to values corresponding to previous works if
!   a user wishes to try to reproduce those works.
! </DESCRIPTION>

! <DATASET NAME="Black carbon emissions">
!   The black carbon emission dataset is that derived in Cooke et al. (1999)
!   The dataset can be obtained from the contact person above.
! </DATASET>

! <DATASET NAME="Organic carbon emissions">
!   The organic carbon emission dataset is that derived in Cooke et al. (1999)
!   The dataset can be obtained from the contact person above.
! </DATASET>

! <INFO>

!   <REFERENCE>
!Cooke, W. F. and J. J. N. Wilson,  A global black carbon aerosol model,
!J. Geophys. Res., 101, 19395-19409, 1996.
! </REFERENCE>
!   <REFERENCE>
!Cooke, W. F., C. Liousse, H. Cachier and J. Feichter,
!Construction of a 1 x 1 fossil fuel emission dataset for carbonaceous
!aerosol and implementation and radiative impact in the ECHAM-4 model,
!J. Geophys. Res., 104, 22137-22162, 1999 </REFERENCE>
!   <REFERENCE>
! Cooke, W.F., V. Ramaswamy and P. Kasibathla,
! A GCM study of the global carbonaceous aerosol distribution.
!J. Geophys. Res., 107, accepted, 2002
! </REFERENCE>
! </INFO>
use types_mod, only : r8
use              fms_mod, only : file_exist,           &
                                 close_file,           &
                                 stdlog,               &
                                 write_version_number
use     time_manager_mod, only : time_type
!use     diag_manager_mod, only : send_data,            &
!                                 register_diag_field,  &
!                                 register_static_field
use   tracer_manager_mod, only : get_tracer_index, &
                                 set_tracer_atts
use    field_manager_mod, only : MODEL_ATMOS
use atmos_tracer_utilities_mod, only : interp_emiss
use        constants_mod, only : PI

implicit none
private
!-----------------------------------------------------------------------
!----- interfaces -------

public  atmos_blackc_sourcesink,   &
        atmos_organic_sourcesink,  &
        atmos_carbon_aerosol_init, &
        atmos_carbon_aerosol_end

!-----------------------------------------------------------------------
!----------- namelist -------------------
!-----------------------------------------------------------------------
!
!  When initializing additional tracers, the user needs to make the
!  following changes.
!
!  Add an integer variable below for each additional tracer. This should
!  be initialized to zero.
!
!  Add id_tracername for each additional tracer. These are used in
!  initializing and outputting the tracer fields.
!
!-----------------------------------------------------------------------

! tracer number for radon
integer :: nbcphobic=0
integer :: nbcphilic=0
integer :: nocphobic=0
integer :: nocphilic=0

!--- identification numbers for  diagnostic fields and axes ----

integer :: id_emissoc, id_emissbc

!--- Arrays to help calculate tracer sources/sinks ---
real(r8), allocatable, dimension(:,:) :: bcsource,ocsource

character(len=6), parameter :: module_name = 'tracer'

logical :: module_is_initialized = .FALSE.
logical :: used

!---- version number -----
character(len=128) :: version = '$Revision$'
character(len=128) :: tagname = '$Id$'
!-----------------------------------------------------------------------

contains

!#######################################################################

!<SUBROUTINE NAME ="atmos_blackc_sourcesink">
!<OVERVIEW>
!  A subroutine to calculate the source and sinks of black carbon aerosol.
!</OVERVIEW>
!
!<DESCRIPTION>
!
! This routine calculates the source and sink terms for black carbon.
! Simply put, the hydrophobic aerosol has sources from emissions and
! sinks from dry deposition and transformation into hydrophilic aerosol.
! The hydrophilic aerosol also has emission sources and has sinks of wet
! and dry deposition.
!
! The following schematic shows how the black carbon scheme
! is implemented.

!

! +------------+  Trans-   +------------+
! |  Hydro-    | formation |  Hydro-    |
! |  phobic    |           |  philic    |
! |  black     |---------->|  black     |
! |  carbon    |           |  carbon    |
! |            |           |            |
! +------------+           +------------+
!    ^      |                ^    |   |
!    |      |                |    |   |
!    |      =                |    =   =
!  Source  Dry            Source Dry Wet
!          Dep.                  Dep Dep
!
!

! The transformation time used here is 1 day, which corresponds to an
! e-folding time of 1.44 days. This can be varied as necessary.

!</DESCRIPTION>
!<TEMPLATE>
!call atmos_blackc_sourcesink (lon, lat, land, pwt, &
!                         black_cphob, black_cphob_dt,  &
!                         black_cphil, black_cphil_dt,  &
!                         Time, is, ie, js, je, kbot)
!</TEMPLATE>
!   <IN NAME="lon" TYPE="real(r8)" DIM="(:,:)">
!     Longitude of the centre of the model gridcells
!   </IN>
!   <IN NAME="lat" TYPE="real(r8)" DIM="(:,:)">
!     Latitude of the centre of the model gridcells
!   </IN>
!   <IN NAME="land" TYPE="real(r8)" DIM="(:,:)">
!     Land/sea mask.
!   </IN>
!   <IN NAME="pwt" TYPE="real(r8)" DIM="(:,:,:)">
!     The pressure weighting array. = dP/grav
!   </IN>
!   <IN NAME="black_cphob" TYPE="real(r8)" DIM="(:,:,:)">
!     The array of the hydrophobic black carbon aerosol mixing ratio
!   </IN>
!   <IN NAME="black_cphil" TYPE="real(r8)" DIM="(:,:,:)">
!     The array of the hydrophilic black carbon aerosol mixing ratio
!   </IN>
!   <IN NAME="Time" TYPE="type(time_type)">
!     Model time.
!   </IN>
!   <IN NAME="is, ie, js, je" TYPE="integer">
!     Local domain boundaries.
!   </IN>
!   <IN NAME="kbot" TYPE="integer, optional" DIM="(:,:)">
!     Integer array describing which model layer intercepts the surface.
!   </IN>

!   <OUT NAME="black_cphob_dt" TYPE="real(r8)" DIM="(:,:,:)">
!     The array of the tendency of the hydrophobic black carbon aerosol mixing ratio.
!   </OUT>
!   <OUT NAME="black_cphil_dt" TYPE="real(r8)" DIM="(:,:,:)">
!     The array of the tendency of the hydrophilic black carbon aerosol mixing ratio.
!   </OUT>

 subroutine atmos_blackc_sourcesink (lon, lat, land, pwt, &
                               black_cphob, black_cphob_dt,  &
                               black_cphil, black_cphil_dt,  &
                               Time, is, ie, js, je, kbot)

!-----------------------------------------------------------------------
   real(r8), intent(in),  dimension(:,:)   :: lon, lat
   real(r8), intent(in),  dimension(:,:)   :: land
   real(r8), intent(in),  dimension(:,:,:) :: pwt, black_cphob,black_cphil
   real(r8), intent(out), dimension(:,:,:) :: black_cphob_dt,black_cphil_dt
type(time_type), intent(in)            :: Time
integer, intent(in)                    :: is, ie, js, je
integer, intent(in),  dimension(:,:), optional :: kbot
!-----------------------------------------------------------------------
   real(r8), dimension(size(black_cphob,1),size(black_cphob,2),size(black_cphob,3)) ::  &
         sourcephob, sinkphob, sourcephil, sinkphil
   real(r8)  dtr
integer  i,j,kb,id,jd,kd,lat1
!-----------------------------------------------------------------------

      id=size(black_cphob,1); jd=size(black_cphob,2); kd=size(black_cphob,3)

      dtr= PI/180.

!----------- compute black carbon source ------------

      sourcephob = 0.0
      sourcephil = 0.0

          do j=1,jd
           sourcephob(:,j,kd)=0.8*bcsource(:,j+js-1)/pwt(:,j,kd)
           sourcephil(:,j,kd)=0.2*bcsource(:,j+js-1)/pwt(:,j,kd) +&
                              8.038e-6*black_cphob(:,j,kd)
          enddo

!------- compute black carbon phobic sink --------------
!
!  BCphob has a half-life time of 1.0days
!   (corresponds to an e-folding time of 1.44 days)
!
!  sink = 1./(86400.*1.44) = 8.023e-6
!

    where (black_cphob(:,:,:) >= 0.0)
       sinkphob(:,:,:) = -8.038e-6*black_cphob(:,:,:)
    elsewhere
       sinkphob(:,:,:) = 0.0
    endwhere

       sinkphil(:,:,:) = 0.0

!------- tendency ------------------

      black_cphob_dt=sourcephob+sinkphob
      black_cphil_dt=sourcephil+sinkphil

!-----------------------------------------------------------------------

 end subroutine atmos_blackc_sourcesink
!</SUBROUTINE >

!#######################################################################

!<SUBROUTINE NAME ="atmos_organic_sourcesink">
!<OVERVIEW>
!  A subroutine to calculate the source and sinks of organic carbon aerosol.
!</OVERVIEW>
!<DESCRIPTION>

! This routine calculates the source and sink terms for organic carbon.
! Simply put, the hydrophobic aerosol has sources from emissions and
! sinks from dry deposition and transformation into hydrophilic aerosol.
! The hydrophilic aerosol also has emission sources and has sinks of wet
! and dry deposition.
!
! The following schematic shows how the organic carbon scheme
! is implemented.

!

! +------------+  Trans-   +------------+
! |  Hydro-    | formation |  Hydro-    |
! |  phobic    |           |  philic    |
! |  organic   |---------->|  organic   |
! |  carbon    |           |  carbon    |
! |            |           |            |
! +------------+           +------------+
!    ^      |                ^    |   |
!    |      |                |    |   |
!    |      =                |    =   =
!  Source  Dry            Source Dry Wet
!          Dep.                  Dep Dep
!

!
! The transformation time used here is 2 days, which corresponds to an
! e-folding time of 2.88 days. This can be varied as necessary.
!
!</DESCRIPTION>
!<TEMPLATE>
!call atmos_organic_sourcesink (lon, lat, land, pwt, organic_carbon, organic_carbon_dt,  &
!                              Time, is, ie, js, je, kbot)
!</TEMPLATE>
!   <IN NAME="lon" TYPE="real(r8)" DIM="(:,:)">
!     Longitude of the centre of the model gridcells
!   </IN>
!   <IN NAME="lat" TYPE="real(r8)" DIM="(:,:)">
!     Latitude of the centre of the model gridcells
!   </IN>
!   <IN NAME="land" TYPE="real(r8)" DIM="(:,:)">
!     Land/sea mask.
!   </IN>
!   <IN NAME="pwt" TYPE="real(r8)" DIM="(:,:,:)">
!     The pressure weighting array. = dP/grav
!   </IN>
!   <IN NAME="organic_carbon" TYPE="real(r8)" DIM="(:,:,:)">
!     The array of the organic carbon aerosol mixing ratio
!   </IN>
!   <IN NAME="Time" TYPE="type(time_type)">
!     Model time.
!   </IN>
!   <IN NAME="is, ie, js, je" TYPE="integer">
!     Local domain boundaries.
!   </IN>
!   <IN NAME="kbot" TYPE="integer, optional" DIM="(:,:)">
!     Integer array describing which model layer intercepts the surface.
!   </IN>

!   <OUT NAME="organic_carbon_dt" TYPE="real(r8)" DIM="(:,:,:)">
!     The array of the tendency of the organic carbon aerosol mixing ratio.
!   </OUT>

 subroutine atmos_organic_sourcesink (lon, lat, land, pwt, organic_carbon, organic_carbon_dt,  &
                              Time, is, ie, js, je, kbot)

!-----------------------------------------------------------------------
   real(r8), intent(in),  dimension(:,:)   :: lon, lat
   real(r8), intent(in),  dimension(:,:)   :: land
   real(r8), intent(in),  dimension(:,:,:) :: pwt, organic_carbon
   real(r8), intent(out), dimension(:,:,:) :: organic_carbon_dt
     type(time_type), intent(in) :: Time
integer, intent(in)                    :: is, ie, js, je
integer, intent(in),  dimension(:,:), optional :: kbot
!-----------------------------------------------------------------------
   real(r8), dimension(size(organic_carbon,1),size(organic_carbon,2),size(organic_carbon,3)) ::  &
         source, sink
   real(r8)  dtr
integer  i,j,kb,id,jd,kd,lat1
!-----------------------------------------------------------------------

      id=size(organic_carbon,1); jd=size(organic_carbon,2); kd=size(organic_carbon,3)

      dtr=PI/180.

!----------- compute organic carbon source ------------

      source = 0.0

      if (present(kbot)) then
          do j=1,jd
          do i=1,id
             kb=kbot(i,j)
             source(i,j,kb)=ocsource(i,j+js-1)/pwt(i,j,kb)
          enddo
          enddo
      else
          do j=1,je-js+1
           source(:,j,kd)= ocsource(:,j+js-1)/pwt(:,j,kd)
          enddo
      endif

!------- compute organic carbon sink --------------
!
!  OCphob has a half-life time of 2.0days
!   (corresponds to an e-folding time of 2.88 days)
!
!  sink = 1./(86400.*2.88) = 4.019e-6
!

    where (organic_carbon(:,:,:) >= 0.0)
         sink(:,:,:) = -4.019e-6*organic_carbon(:,:,:)
    elsewhere
       sink(:,:,:) = 0.0
    endwhere

!------- tendency ------------------

      organic_carbon_dt=source+sink

!-----------------------------------------------------------------------

 end subroutine atmos_organic_sourcesink
!</SUBROUTINE>

!#######################################################################

!<SUBROUTINE NAME ="atmos_carbon_aerosol_init">

!<OVERVIEW>
! Subroutine to initialize the carbon aerosol module.
!</OVERVIEW>
!<DESCRIPTION>
! This subroutine querys the tracer manager to find the indices for the
! various carbonaceous aerosol tracers. It also registers the emission
! fields for diagnostic purposes.
!
!</DESCRIPTION>
!<TEMPLATE>
!call atmos_carbon_aerosol_init (lonb, latb, r, axes, Time, mask)
!</TEMPLATE>
!   <IN NAME="lonb" TYPE="real(r8)" DIM="(:)">
!     The longitudes for the local domain.
!   </IN>
!   <IN NAME="latb" TYPE="real(r8)" DIM="(:)">
!     The latitudes for the local domain.
!   </IN>
!   <INOUT NAME="r" TYPE="real(r8)" DIM="(:,:,:,:)">
!     Tracer fields dimensioned as (nlon,nlat,nlev,ntrace).
!   </INOUT>
!   <IN NAME="mask" TYPE="real(r8), optional" DIM="(:,:,:)">
!      optional mask (0. or 1.) that designates which grid points
!           are above (=1.) or below (=0.) the ground dimensioned as
!           (nlon,nlat,nlev).
!   </IN>
!   <IN NAME="Time" TYPE="type(time_type)">
!     Model time.
!   </IN>
!   <IN NAME="axes" TYPE="integer" DIM="(4)">
!     The axes relating to the tracer array dimensioned as
!      (nlon, nlat, nlev, ntime)
!   </IN>

 subroutine atmos_carbon_aerosol_init (lonb, latb, r, axes, Time, mask)

!-----------------------------------------------------------------------
!
!   r    = tracer fields dimensioned as (nlon,nlat,nlev,ntrace)
!   mask = optional mask (0. or 1.) that designates which grid points
!          are above (=1.) or below (=0.) the ground dimensioned as
!          (nlon,nlat,nlev).
!
!-----------------------------------------------------------------------
real(r8), dimension(:),    intent(in) :: lonb, latb
real(r8),            intent(inout), dimension(:,:,:,:) :: r
integer        , intent(in)                        :: axes(4)
type(time_type), intent(in)                        :: Time
real(r8),            intent(in),    dimension(:,:,:), optional :: mask

integer :: n

   if (module_is_initialized) return

!----- set initial value of carbon ------------

   n = get_tracer_index(MODEL_ATMOS,'bcphob')
   if (n>0) then
      nbcphobic = n
      call set_tracer_atts(MODEL_ATMOS,'bcphob','hphobic_bc','g/g')
      if (nbcphobic > 0 ) write (*,30) 'Hydrophobic BC',nbcphobic
      if (nbcphobic > 0 ) write (stdlog(),30) 'Hydrophobic BC',nbcphobic
   endif

   n = get_tracer_index(MODEL_ATMOS,'bcphil')
   if (n>0) then
      nbcphilic=n
      call set_tracer_atts(MODEL_ATMOS,'bcphil','hphilic_bc','g/g')
      if (nbcphilic > 0 ) write (*,30) 'Hydrophilic BC',nbcphilic
      if (nbcphilic > 0 ) write (stdlog(),30) 'Hydrophilic BC',nbcphilic
   endif

   n = get_tracer_index(MODEL_ATMOS,'ocphob')
   if (n>0) then
      nocphobic=n
      call set_tracer_atts(MODEL_ATMOS,'ocphob','hphobic_oc','g/g')
      if (nocphobic > 0 ) write (*,30) 'Hydrophobic OC',nocphobic
      if (nocphobic > 0 ) write (stdlog(),30) 'Hydrophobic OC',nocphobic
   endif

   n = get_tracer_index(MODEL_ATMOS,'ocphil')
   if (n>0) then
      nocphilic=n
      call set_tracer_atts(MODEL_ATMOS,'ocphil','hphilic_oc','g/g')
      if (nocphilic > 0 ) write (*,30) 'Hydrophilic OC',nocphilic
      if (nocphilic > 0 ) write (stdlog(),30) 'Hydrophilic OC',nocphilic
   endif

  30        format (A,' was initialized as tracer number ',i2)
      !Read in emission files
!
!   id_emissbc = register_static_field ( 'tracers',                    &
!                     'bcemiss', axes(1:2),       &
!                     'bcemiss', 'g/m2/s')
!   id_emissoc = register_static_field ( 'tracers',                    &
!                     'ocemiss', axes(1:2),       &
!                     'ocemiss', 'g/m2/s')
!
   allocate (bcsource(size(lonb)-1,size(latb)-1))
   allocate (ocsource(size(lonb)-1,size(latb)-1))
   call tracer_input(lonb, latb, Time)

   call write_version_number (version, tagname)
   module_is_initialized = .TRUE.

!-----------------------------------------------------------------------

end subroutine atmos_carbon_aerosol_init
!</SUBROUTINE>

!<SUBROUTINE NAME ="atmos_carbon_aerosol_end">
!<OVERVIEW>
!  The destructor routine for the carbon aerosol module.
!</OVERVIEW>
! <DESCRIPTION>
! This subroutine writes the version name to logfile and exits.
! </DESCRIPTION>
!<TEMPLATE>
!call atmos_carbon_aerosol_end
!</TEMPLATE>
subroutine atmos_carbon_aerosol_end

   module_is_initialized = .FALSE.

 end subroutine atmos_carbon_aerosol_end
!</SUBROUTINE>

!#######################################################################
 subroutine tracer_input(lonb, latb, Time)
real(r8), dimension(:),    intent(in) :: lonb, latb
type(time_type),intent(in) :: Time

integer      :: i, j, unit, io
real(r8)         :: emiss
real(r8)         :: dtr, deg_90, deg_180, deg3p6, deg3!, modxdeg, modydeg
real(r8)         :: ZCARBONSEASON(12)
real(r8)         :: bcsource1(100,60)
logical :: opened
!
! This is the Rotty seaonality for fossil fuel emissions of sulfate.
!
      DATA ZCARBONSEASON/1.146,1.139,1.081,0.995,0.916,0.920,0.910, &
                         0.907,0.934,0.962,1.019,1.072/
!
      dtr= PI/180.
      deg_90= -90.*dtr; deg_180= -180.*dtr
      ! -90 and -180 degrees are the southwest boundaries of the
      ! emission  field you are reading in.
      deg3p6 = 3.6*dtr; deg3 = 3.*dtr;
      ! 3.6 degrees longitude and 3 degree latitude is the resolution
      ! of the r30 emission data that I used in SKYHI.

! initialise the BC phobic and philic
! read in the emission sources here.

do unit = 30,100
INQUIRE(unit=unit, opened= opened)
if (.NOT. opened) exit
enddo
         open (unit,file='INPUT/r30.bc.ann', form='formatted', action='read')
        do io= 1,6000
           read  (unit, FMT=1968, end=11) i,j,emiss
           bcsource1(i,j)=emiss
        enddo
  11    call close_file (unit)
1968  FORMAT(2I3,e11.4)
1969  FORMAT(2I3,f11.3)
! Interpolate the R30 emission field to the resolution of the model.
        call interp_emiss ( bcsource1, 0.0_r8, deg_90, deg3p6, deg3, &
                     bcsource)

 write(*,*) 'Reading OC emissions'
!Now let's do the OC
!
         bcsource1 = 0.0E+00
do unit = 30,100
INQUIRE(unit=unit, opened= opened)
if (.NOT. opened) exit
enddo
         open (unit,file='INPUT/r30.oc.ann', form='formatted', action='read')
         do io= 1,6000
            read  (unit, FMT=1968, end=13) i,j,emiss
            bcsource1(i,j)=emiss
         enddo
  13     call close_file (unit)

! Interpolate the R30 emission field to the resolution of the model.
         call interp_emiss ( bcsource1, 0.0_r8, deg_90, deg3p6, deg3, &
                              ocsource)

! Send the emission data to the diag_manager for output.
!         if (id_emissbc > 0 ) &
!           used = send_data ( id_emissbc, bcsource, Time )
!         if (id_emissoc > 0 ) &
!           used = send_data ( id_emissoc, ocsource, Time )

end subroutine tracer_input

end module atmos_carbon_aerosol_mod

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
!
! Sensor_Planck_Utility module
!
! Container module for all the Planck and Sensor_Planck modules and dependencies.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 16-Feb-2012
!                       paul.vandelst@noaa.gov
!

MODULE Sensor_Planck_Utility
  ! The Utility modules
  USE Type_Kinds
  USE File_Utility
  USE Message_Handler
  USE Fundamental_Constants
  USE Compare_Float_Numbers
  ! The SpcCoeff modules
  USE SpcCoeff_Define
  USE SpcCoeff_Binary_IO
  ! The Planck modules
  USE Planck_Functions
  USE Sensor_Planck_Functions
  PUBLIC
END MODULE Sensor_Planck_Utility

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
! Like vector_subscript_1.f90, but check subscripts in multi-dimensional
! arrays.
! { dg-do run }
program main
  implicit none
  integer, parameter :: n = 5
  integer :: i1, i2, i3
  integer, dimension (n, n, n) :: a, b
  integer, dimension (n) :: idx, id

  idx = (/ 3, 1, 5, 2, 4 /)
  id = (/ (i1, i1 = 1, n) /)
  forall (i1 = 1:n, i2 = 1:n, i3 = 1:n)
    b (i1, i2, i3) = i1 + i2 * 10 + i3 * 100
  end forall

  i1 = 5
  a (foo (i1), 1, :) = b (2, :, foo (i1))
  do i1 = 1, 5
    do i2 = 1, 5
      if (a (idx (i1), 1, i2) .ne. b (2, i1, idx (i2))) call abort
    end do
  end do
  a = 0

  a (1, idx (1:4), 2:4) = b (2:5, idx (3:5), 2)
  do i1 = 1, 4
    do i2 = 1, 3
      if (a (1, idx (i1), 1 + i2) .ne. b (1 + i1, idx (i2 + 2), 2)) call abort
    end do
  end do
  a = 0
contains
  function foo (n)
    integer :: n
    integer, dimension (n) :: foo
    foo = idx (1:n)
  end function foo
end program main

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text: #redirect wiki:Sacramento:Red Sea Food Market and Halal Meat

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*
* $Id$
*

      subroutine integrate_kbppv3d(version,rlocal,
     >                            nrho,drho,lmax,locp,zv,
     >                            vp,wp,rho,f,cs,sn,
     >                            nfft3d,lmmax,
     >                            G,vl,vnl,
     >                            n_prj,l_prj,m_prj,b_prj,vnlnrm,
     >                            semicore,rho_sc_r,rho_sc_k,
     >                            ierr)
      implicit none
      integer          version
      double precision rlocal
      integer          nrho
      double precision drho
      integer          lmax
      integer          locp
      double precision zv
      double precision vp(nrho,0:lmax)
      double precision wp(nrho,0:lmax)
      double precision rho(nrho)
      double precision f(nrho)
      double precision cs(nrho)
      double precision sn(nrho)

      integer nfft3d,lmmax
      double precision G(nfft3d,3)
      double precision vl(nfft3d)
      double precision vnl(nfft3d,lmmax)
      integer          n_prj(lmmax),l_prj(lmmax),m_prj(lmmax)
      integer          b_prj(lmmax)
      double precision vnlnrm(0:lmax)

      logical semicore
      double precision rho_sc_r(nrho,2)
      double precision rho_sc_k(nfft3d,4)

      integer ierr

      integer np,taskid,MASTER
      integer np_i,np_j,taskid_i,taskid_j,countj
      parameter (MASTER=0)

*     *** local variables ****
      logical fast_erf
      integer lcount
      integer k1,k2,k3,i,l,pzero,zero
      double precision pi,twopi,forpi
      double precision p0,p1,p2,p3,p
      double precision gx,gy,gz,a,q,d
      double precision tollz
      parameter(tollz=1d-16)

*     **** Error function parameters ****
      real*8 yerf,xerf
      real*8 c1,c2,c3,c4,c5,c6
      parameter (c1=0.07052307840d0,c2=0.04228201230d0)
      parameter (c3=0.00927052720d0)
      parameter (c4=0.00015201430d0,c5=0.00027656720d0)
      parameter (c6=0.00004306380d0)

*     **** external functions ****
      logical          control_fast_erf
      external         control_fast_erf
      double precision dsum,simp,util_erf
      external         dsum,simp,util_erf

      call Parallel2d_np_i(np_i)
      call Parallel2d_np_j(np_j)
      call Parallel2d_taskid_i(taskid_i)
      call Parallel2d_taskid_j(taskid_j)

      fast_erf = control_fast_erf()

      pi=4.0d0*datan(1.0d0)
      twopi=2.0d0*pi
      forpi=4.0d0*pi

      IF(LMMAX.GT.16) THEN
        IERR=1
        RETURN
      ENDIF
      IF((NRHO/2)*2.EQ.NRHO) THEN
        IERR=2
        RETURN
      ENDIF

      P0=DSQRT(FORPI)
      P1=DSQRT(3.0d0*FORPI)
      P2=DSQRT(15.0d0*FORPI)
      P3=DSQRT(105.0d0*FORPI)

*::::::::::::::::::  Define non-local pseudopotential  ::::::::::::::::
      do l=0,lmax
        if (l.ne.locp) then
          do I=1,nrho
            vp(i,l)=vp(i,l)-vp(i,locp)
          end do
        end if
      end do

*:::::::::::::::::::::  Normarization constants  ::::::::::::::::::::::
      lcount = 0
      do l=0,lmax
        if (l.ne.locp) then
          do i=1,nrho
            f(i)=vp(i,l)*wp(i,l)**2
          end do
          a=simp(nrho,f,drho)
          vnlnrm(l) = (1.0d0/a)
        else
          vnlnrm(l) = 0.0d0
        end if
      end do

*======================  Fourier transformation  ======================
      call dcopy(nfft3d,0.0d0,0,vl,1)
      call dcopy(lmmax*nfft3d,0.0d0,0,vnl,1)
      call dcopy(4*nfft3d,0.0d0,0,rho_sc_k,1)

*     ***** find the G==0 point in the lattice *****
      call D3dB_ijktoindexp(1,1,1,1,zero,pzero)

      countj = -1
      DO 700 k1=1,nfft3d

        countj = mod(countj+1,np_j)

        if (countj.ne.taskid_j) go to 700
        if ((pzero.eq.taskid_i).and.(k1.eq.zero)) go to 700

        Q=DSQRT(G(k1,1)**2
     >         +G(k1,2)**2
     >         +G(k1,3)**2)
        if (abs(q).lt.tollz) go to 700

        GX=G(k1,1)/Q
        GY=G(k1,2)/Q
        GZ=G(k1,3)/Q
        DO I=1,NRHO
          CS(I)=DCOS(Q*RHO(I))
          SN(I)=DSIN(Q*RHO(I))
        END DO

        lcount = lmmax+1
        GO TO (500,400,300,200), LMAX+1

*::::::::::::::::::::::::::::::  f-wave  ::::::::::::::::::::::::::::::
  200   CONTINUE
        if (locp.ne.3) then
           F(1)=0.0d0
           do I=2,NRHO
             A=SN(I)/(Q*RHO(I))
             A=15.0d0*(A-CS(I))/(Q*RHO(I))**2 - 6*A + CS(I)
             F(I)=A*WP(I,3)*VP(I,3)
           end do
           D=P3*SIMP(NRHO,F,DRHO)/Q

           lcount = lcount-1
           vnl(k1,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY)
     >                          /dsqrt(24.0d0)
           lcount = lcount-1
           vnl(k1,lcount)=D*GX*GY*GZ
           lcount = lcount-1
           vnl(k1,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0)
     >                          /dsqrt(40.0d0)
           lcount = lcount-1
           vnl(k1,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0)
     >                          /dsqrt(60.0d0)
           lcount = lcount-1
           vnl(k1,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0)
     >                          /dsqrt(40.0d0)
           lcount = lcount-1
           vnl(k1,lcount)=D*GZ*(GX*GX - GY*GY)
     >                          /2.0d0
           lcount = lcount-1
           vnl(k1,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ))
     >                          /dsqrt(24.0d0)
c           lcount = lcount-1
c           vnl(k1,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ))
c     >                          /dsqrt(24.0d0)
c           lcount = lcount-1
c           vnl(k1,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY)
c     >                          /dsqrt(24.0d0)
c           lcount = lcount-1
c           vnl(k1,lcount)=D*GZ*(GX*GX - GY*GY)
c     >                          /2.0d0
c           lcount = lcount-1
c           vnl(k1,lcount)=D*GX*GY*GZ
c           lcount = lcount-1
c           vnl(k1,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0)
c     >                          /dsqrt(40.0d0)
c           lcount = lcount-1
c           vnl(k1,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0)
c     >                          /dsqrt(40.0d0)
c           lcount = lcount-1
c           vnl(k1,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0)
c     >                          /dsqrt(60.0d0)
        end if

*::::::::::::::::::::::::::::::  d-wave  ::::::::::::::::::::::::::::::
  300   CONTINUE
        if (locp.ne.2) then
          F(1)=0.0d0
          DO I=2,NRHO
            A=3.0d0*(SN(I)/(Q*RHO(I))-CS(I))/(Q*RHO(I))-SN(I)
            F(I)=A*WP(I,2)*VP(I,2)
          END DO
          D=P2*SIMP(NRHO,F,DRHO)/Q

          lcount = lcount-1
          vnl(k1,lcount)=D*GX*GY
          lcount = lcount-1
          vnl(k1,lcount)=D*GY*GZ
          lcount = lcount-1
          vnl(k1,lcount)=D*(3.0d0*GZ*GZ-1.0d0)
     >                          /(2.0d0*dsqrt(3.0d0))
          lcount = lcount-1
          vnl(k1,lcount)=D*GZ*GX
          lcount = lcount-1
          vnl(k1,lcount)=D*(GX*GX-GY*GY)/(2.0d0)

c          lcount = lcount-1
c          vnl(k1,lcount)=D*(3.0d0*GZ*GZ-1.0d0)
c     >                          /(2.0d0*dsqrt(3.0d0))
c          lcount = lcount-1
c          vnl(k1,lcount)=D*GX*GY
c          lcount = lcount-1
c          vnl(k1,lcount)=D*GY*GZ
c          lcount = lcount-1
c          vnl(k1,lcount)=D*GZ*GX
c          lcount = lcount-1
c          vnl(k1,lcount)=D*(GX*GX-GY*GY)/(2.0d0)
        end if

*::::::::::::::::::::::::::::::  p-wave  ::::::::::::::::::::::::::::::
  400   CONTINUE
        if (locp.ne.1) then
           F(1)=0.0d0
           DO I=2,NRHO
             F(I)=(SN(I)/(Q*RHO(I))-CS(I))*WP(I,1)*VP(I,1)
           END DO
           P=P1*SIMP(NRHO,F,DRHO)/Q
           lcount = lcount-1
           vnl(k1,lcount)=P*GY
           lcount = lcount-1
           vnl(k1,lcount)=P*GZ
           lcount = lcount-1
           vnl(k1,lcount)=P*GX
        end if

*::::::::::::::::::::::::::::::  s-wave  :::::::::::::::::::::::::::::::
  500   CONTINUE
        if (locp.ne.0) then
          DO I=1,NRHO
            F(I)=SN(I)*WP(I,0)*VP(I,0)
          END DO
          lcount = lcount-1
          vnl(k1,lcount)=P0*SIMP(NRHO,F,DRHO)/Q
        end if

*::::::::::::::::::::::::::::::  local  :::::::::::::::::::::::::::::::
  600   CONTINUE

        if (version.eq.3) then
        DO  I=1,NRHO
          F(I)=RHO(I)*VP(I,locp)*SN(I)
        END DO
        vl(k1)=SIMP(NRHO,F,DRHO)*FORPI/Q-ZV*FORPI*CS(NRHO)/(Q*Q)
        end if

        if (version.eq.4) then
        if (fast_erf) then
           do I=1,NRHO
             xerf=RHO(I)/rlocal
             yerf = (1.0d0
     >            + xerf*(c1 + xerf*(c2
     >            + xerf*(c3 + xerf*(c4
     >            + xerf*(c5 + xerf*c6))))))**4
             yerf = (1.0d0 - 1.0d0/yerf**4)
             F(I)=(RHO(I)*VP(I,locp)+ZV*yerf)*SN(I)
           end do
        else
           do I=1,NRHO
             xerf=RHO(I)/rlocal
             yerf = util_erf(xerf)
             F(I)=(RHO(I)*VP(I,locp)+ZV*yerf)*SN(I)
           end do
        end if
        vl(k1)=SIMP(NRHO,F,DRHO)*FORPI/Q
        end if

*::::::::::::::::::::: semicore density :::::::::::::::::::::::::::::::
        if (semicore) then
           do i=1,nrho
              f(i) = rho(i)*dsqrt(rho_sc_r(i,1))*sn(i)
           end do
           rho_sc_k(k1,1) = SIMP(nrho,f,drho)*forpi/Q

           do i=1,nrho
             f(i)=(sn(i)/(Q*rho(i))-cs(i))*rho_sc_r(i,2)*rho(i)
           end do
           P = SIMP(nrho,f,drho)*forpi/Q
           rho_sc_k(k1,2)=P*GX
           rho_sc_k(k1,3)=P*GY
           rho_sc_k(k1,4)=P*GZ

        end if

  700 CONTINUE
      call D1dB_Vector_SumAll(4*nfft3d,rho_sc_k)
      call D1dB_Vector_SumAll(nfft3d,vl)
      call D1dB_Vector_Sumall(lmmax*nfft3d,vnl)

*:::::::::::::::::::::::::::::::  G=0  ::::::::::::::::::::::::::::::::
      if (pzero.eq.taskid_i) then

         if (version.eq.3) then
         DO I=1,NRHO
           F(I)=VP(I,locp)*RHO(I)**2
         END DO
         vl(zero)=FORPI*SIMP(NRHO,F,DRHO)+TWOPI*ZV*RHO(NRHO)**2
         end if

         if (version.eq.4) then
         if (fast_erf) then
            do I=1,NRHO
              xerf=RHO(I)/rlocal
              yerf = (1.0d0
     >             + xerf*(c1 + xerf*(c2
     >             + xerf*(c3 + xerf*(c4
     >             + xerf*(c5 + xerf*c6))))))**4
              yerf = (1.0d0 - 1.0d0/yerf**4)
              F(I)=(VP(I,locp)*RHO(I)+ZV*yerf)*RHO(I)
            end do
         else
            do I=1,NRHO
              xerf=RHO(I)/rlocal
              yerf = util_erf(xerf)
              F(I)=(VP(I,locp)*RHO(I)+ZV*yerf)*RHO(I)
            end do
         end if
         vl(zero)=FORPI*SIMP(NRHO,F,DRHO)
         end if

*        **** semicore density ****
         if (semicore) then
            do i=1,nrho
               f(i) = dsqrt(rho_sc_r(i,1))*rho(i)**2
            end do
            rho_sc_k(zero,1) = forpi*SIMP(nrho,f,drho)
            rho_sc_k(zero,2) = 0.0d0
            rho_sc_k(zero,3) = 0.0d0
            rho_sc_k(zero,4) = 0.0d0
         end if

         do l=1,lmmax
           vnl(zero,l)=0.0d0
         end do
*        *** only j0 is non-zero at zero ****
         if (locp.ne.0) then
            DO  I=1,NRHO
              F(I)=RHO(I)*WP(I,0)*VP(I,0)
            END DO
            vnl(zero,1)=P0*SIMP(NRHO,F,DRHO)
         end if

      end if

*     ********************************
*     **** define n_prj and l_prj ****
*     ********************************
      lcount = lmmax+1
      GO TO (950,940,930,920), lmax+1

        !::::::  f-wave  :::::::
  920   CONTINUE
        if (locp.ne.3) then
          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 3
          m_prj(lcount) = -3
          b_prj(lcount) = 4

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 3
          m_prj(lcount) = -2
          b_prj(lcount) = 4

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 3
          m_prj(lcount) = -1
          b_prj(lcount) = 4

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 3
          m_prj(lcount) = 0
          b_prj(lcount) = 4

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 3
          m_prj(lcount) = 1
          b_prj(lcount) = 4

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 3
          m_prj(lcount) = 2
          b_prj(lcount) = 4

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 3
          m_prj(lcount) = 3
          b_prj(lcount) = 4
        end if

        !::::  d-wave  ::::
  930   CONTINUE
        if (locp.ne.2) then
          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 2
          m_prj(lcount) = -2
          b_prj(lcount) = 3

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 2
          m_prj(lcount) = -1
          b_prj(lcount) = 3

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 2
          m_prj(lcount) = 0
          b_prj(lcount) = 3

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 2
          m_prj(lcount) = 1
          b_prj(lcount) = 3

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 2
          m_prj(lcount) = 2
          b_prj(lcount) = 3
        end if

        !::::  p-wave  ::::
  940   CONTINUE
        if (locp.ne.1) then
          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 1
          m_prj(lcount) = -1
          b_prj(lcount) = 2

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 1
          m_prj(lcount) = 0
          b_prj(lcount) = 2

          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 1
          m_prj(lcount) = 1
          b_prj(lcount) = 2
        end if

        !::::  s-wave  ::::
  950   CONTINUE
        if (locp.ne.0) then
          lcount = lcount-1
          n_prj(lcount) = 1
          l_prj(lcount) = 0
          m_prj(lcount) = 0
          b_prj(lcount) = 1
        end if

      IERR=0
      RETURN
      END

max_stars_repo_name: masamichitakagi/nwchem

max_issues_repo_name: masamichitakagi/nwchem

max_forks_repo_name: masamichitakagi/nwchem


text:
Lab mice are much like regular mice, but are used as tools in scientific study and experiementation.  The stereotypical lab mouse is an albino, with white fur and skin and red eyes.  Some have been specially bred or genetically engineered to exhibit certain traits which are then studied.

There is considerable controversy over whether this practice is ethical, but it goes on nonetheless.

Lab mice at UC Davis are occasionally sent to other universities for study there.  Few transport companies will accept live animals; DHL is one of the few that does.  One of the worst parts of being an employee in the Central Stores and Receiving department at UC Davis is periodically having to listen to 100+ tiny lab mice skitter and scratch around in their cardboard shipping boxes until the transport folk arrive  like nails on a chalkboard, only a thousand times worse.

It should be noted that lab mice (and their PetCo equivelent for sale as snake food) tend to make poor pets for several reasons.  They tend to be bred to purpose (and inbred), and as a result display neurotic behaviors.  They also get along poorly with mice that have a better breeding heritage (generally referred to as fancy mice).


text:
C--------------------------------------------------------------------------
        subroutine ftpknj(ounit,keywrd,nstart,nkey,intval,comm,
     &                    status)

C       write an array of integer values to header records
C
C       ounit   i  fortran output unit number
C       keywrd  c  keyword name    ( 8 characters, cols.  1- 8)
C       nstart  i  starting sequence number (usually 1)
C       nkey    i  number of keywords to write
C       intval  i  array of keyword values
C       comm    c  array of keyword comments (47 characters, cols. 34-80)
C       OUTPUT PARAMETERS:
C       status  i  output error status (0 = ok)
C
C       written by Wm Pence, HEASARC/GSFC, June 1991

        character*(*) keywrd,comm(*)
        integer nstart,nkey,ounit,status,intval(*),i,j
        character keynam*8,comm1*48
        logical repeat

        if (status .gt. 0)return

C       check if the first comment string is to be repeated for all keywords
C       (if the last non-blank character is '&', then it is to be repeated)
        call ftcrep(comm(1),comm1,repeat)

        j=nstart
        do 10 i=1,nkey
C               construct keyword name:
                call ftkeyn(keywrd,j,keynam,status)

C               write the keyword record
                if (repeat)then
                   call ftpkyj(ounit,keynam,intval(i),comm1,status)
                else
                   call ftpkyj(ounit,keynam,intval(i),comm(i),status)
                end if
                if (status .gt. 0)return
                j=j+1
10      continue
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
      subroutine lid_cistern(sb,j,k,lid_prec)

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    Simulate cistern processes

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~ ~ ~ ~ ~
!!    name             |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    sb               |none          |Subbasin number
!!    j                |none          |HRU number
!!    k                |none          |Subdaily time index
!!    lid_prec         |mm            |Precipitation depth a LID receives in a simulation time interval
!!    idt              |minutes       |Simulation time interval for sub-daily modeling
!!    ihru             |none          |HRU number
!!    nstep            |none          |Number of time intervals for a day
!!    urblu(:)         |none          |Urban land type identification number from urban.dat
!!    lid_vol          |m3            |Volume of a cistern storage
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    name             |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    lid_qsurf        |mm H2O        |Depth of runoff generated on a LID in a given time interval
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name             |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    jj               |none          |Urban land type identification number from
!!                                    |urban.dat
!!    dt               |hour          |Time interval in hours
!!    lid_cumr         |mm H2O        |Cumulative amount of rainfall a LID receives in a time interval
!!    Cumulative amount of rainfall a LID receives in a time interval        |mm H2O        |Cumulative amount of excess rainfall at a time step in a day
!!    lid_exinc        |mm H2O        |Amount of excess rainfall
!!    lid_str          |m3            |Volume of stored water in the cistern
!!    lid_vbypass      |m3            |Volume of flow bypassing the cistern
!!    lid_bypass       |mm            |Depth of flow bypassing the cistern
!!    lid_irr          |m3            |Volume of water taken from the cistern for irrigation
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    None

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

      use parm
      implicit none

      integer :: jj,sb,j,k,temp_count
      real*8 :: lid_str,lid_vbypass,lid_bypass,lid_irr,lid_vol,
     & lid_cumirr
      real :: lid_prec

      jj = urblu(j)

      lid_vol = cs_vol(sb,jj)         ! commented out for testing
      if (lid_vol <= 0) then
        lid_vol = (cs_rdepth(sb,jj) / 1000.) *
     &  (lid_farea(j,3) * fcimp(urblu(j)) * hru_ha(j) * 10000.)         ! for testing, assuming the storage can handle excess rainfall of 2.5 mm
      end if

      lid_str = lid_str_last(j,3)
      lid_cumirr = lid_cumirr_last(j,3)
      lid_bypass = 0.

      if (hrnopcp(sb,k-1) > 96) then ! four days
          lid_irr = 0.3 * lid_vol / nstep     ! assumming 30% of water storage of a cistern in a day
      else
          lid_irr = 0.
      end if

      if (cs_grcon(sb,jj)==0) then
        lid_str = lid_str_last(j,3) + (lid_prec / 1000.) *
     &  (lid_farea(j,3) * fcimp(urblu(j)) * hru_ha(j) * 10000.)           ! m3
      else
        lid_str = lid_str_last(j,3) + (lid_prec / 1000.) *
     &  (hru_ha(j) * 10000.)                                              ! m3
      end if

      if (lid_str > lid_vol) then
          lid_vbypass = lid_str - lid_vol ! assuming water stored in a cistern is used
          lid_str = lid_vol
          lid_bypass = lid_vbypass /
     & (hru_ha(j) * 10000.) * 1000.  ! mm
!     & (lid_farea(j,3) * fcimp(urblu(j)) * hru_ha(j) * 10000.) * 1000.  ! mm
      else
          lid_vbypass = 0.
      end if

      lid_qsurf(j,3) = lid_bypass

!! begin temporary
!      if (k == nstep+1) then
!     if (sb == 1) then
!        if (jj == 14) then
!        if (j == 1) then
!          write (3333333,'(5f12.4)') lid_prec,lid_str,lid_irr,
!     & lid_bypass,lid_qsurf(j,3)
!        end if
!      end if
!      end if
!! end temporary

      lid_str = lid_str - lid_irr
      if (lid_str < 0) then
        lid_str = 0.
        lid_irr = lid_str
      end if

      lid_cumirr = lid_cumirr + lid_irr

      if (k == nstep+1) then
          lid_sw_add(j,3) = lid_cumirr / (hru_ha(j) * 10000.) * 1000.
          lid_cumirr = 0.
      end if

      lid_cumirr_last(j,3) = lid_cumirr
      lid_str_last(j,3) = lid_str

      return
      end subroutine