!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
*** S/P EXTDIAG

      subroutine extdiag (d,f,v,dsiz,fsiz,vsiz, 1,48
     $                    trnch,icpu,ni,nk)
*
#include "impnone.cdk"
#include "phy_macros_f.h"
*
      integer dsiz,fsiz,vsiz,icpu,trnch,ni,nk
      real d(dsiz),f(fsiz), v(vsiz)
*
*
*Author
*          B. Bilodeau Feb 2003 - from serdyn5 and phyexe1
*
*
*Revisions
*001      J.P. Toviessi (Aug 2003) - IBM conversion
*               - calls to vssqrt routine (from massvp4 library)
*               - calls to exponen4 (to calculate power function '**')
*               - calls to optimized routine MFOHR,MFOHRA
*
*002      B. Bilodeau (Sep 2003) - exponen4 replaced by vspown1
*003      A. Glazer (September 2003) - Calculation of TD
*
*
*Object
*          to calculate averages and accumulators 
*          of tendencies and diagnostics
*
*Arguments
*
*          - Input/Output -
* F        permanent bus
*
*          - Input -
* D        dynamic bus
* V        volatile (output) bus
*
*          - Input -
* DSIZ     dimension of d
* FSIZ     dimension of f
* VSIZ     dimension of v
* TRNCH    slice number
* ICPU     task number
* N        horizontal running length
* NK       vertical dimension
*
*
*IMPLICITES
*
#include "options.cdk"
#include "sercmdk.cdk"
#include "phybus.cdk"
*
*
*MODULES
*
      EXTERNAL SERXST
      EXTERNAL MVZNXST,MZONXST
*
*
*
**
*
      INTEGER I,K,IK,IERGET
      REAL HEURSER
*
      REAL tmp1
      INTEGER MODP
************************************************************************
*     AUTOMATIC ARRAYS
************************************************************************
*
*
      AUTOMATIC (  P      , REAL , (NI,NK  ) )
      AUTOMATIC (  WORK2D1, REAL , (NI,NK  ) )
      AUTOMATIC (  WORK2D2, REAL , (NI,NK  ) )
      AUTOMATIC (  WORK1D , REAL , (NI     ) )
*

************************************************************************
*
*
#include "consphy.cdk"
*
*     fonction-formule pour faciliter le calcul des indices
      ik(i,k) = (k-1)*ni + i -1


*
      call serget ( 'HEURE' , heurser , 1 , ierget  )
*
*
*     EXTRACTION OF LATITUDE, LONGITUDE AND ARRAY OF ONES
*     (TO VALIDATE ZONAL DIAGNOTICS)
      do i=1,ni
         work1d(i) = 1.0
      end do
*
      call mzonxst(f(dlat) , 'LA', trnch,ni,heurser,1.,  -1, icpu)
      call mzonxst(f(dlon) , 'LO', trnch,ni,heurser,1.,  -1, icpu)
      call mzonxst(work1d  , 'UN', trnch,ni,heurser,1.,  -1, icpu)
*
*
*     EXTRACTION OF ACCUMULATED FLUXES
*
      call serxst  ( f(eiaf) , 'AR', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(eiaf) , 'AR', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(evaf) , 'AU', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(evaf) , 'AU', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(fcaf) , 'AH', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(fcaf) , 'AH', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(fiaf) , 'AD', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(fiaf) , 'AD', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(fqaf) , 'AW', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(fqaf) , 'AW', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(suaf) , 'S7', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(suaf) , 'S7', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(svaf) , 'S8', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(svaf) , 'S8', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(fsaf) , 'AS', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(fsaf) , 'AS', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(fvaf) , 'AV', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(fvaf) , 'AV', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(ivaf) , 'AB', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(ivaf) , 'AB', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(ntaf) , 'NF', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(ntaf) , 'NF', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(siaf) , 'AI', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(siaf) , 'AI', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(flaf) , 'AG', trnch, ni, 0.,      1.,  -1      )
      call mzonxst ( f(flaf) , 'AG', trnch, ni, heurser, 1.,  -1, icpu)
      call serxst  ( f(flusolaf),'AF',trnch,ni, 0.    ,  1.,  -1      )
      call mzonxst ( f(flusolaf),'AF',trnch,ni, heurser, 1.,  -1, icpu)
*
*
*     EXTRACTION OF TENDENCIES DUE TO CONVECTION
*
      call serxst  ( v( tcond), 'TA', trnch, ni, 0., 1.,        -1)
      call mzonxst ( v( tcond), 'TA', trnch, ni, heurser, 
     $               d(pmoins),-2, icpu)
      call serxst  ( v(hucond), 'QA', trnch, ni, 0.,      1.,   -1)
      call mzonxst ( v(hucond), 'QA', trnch, ni, heurser, 
     $               d(pmoins),-2, icpu)
      call serxst  ( v(qccond), 'QT', trnch, ni, 0.,      1.,   -1)
      call mzonxst ( v(qccond), 'QT', trnch, ni, heurser, 
     $               d(pmoins),-2, icpu)
*
*     EXTRACTION OF PRECIPITATION RATES
*
      call serxst (f(tlc)  , 'P1', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(f(tlc)  , 'P1', trnch,ni,heurser,1.,  -1, icpu)
      call serxst (f(tls)  , 'P2', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(f(tls)  , 'P2', trnch,ni,heurser,1.,  -1, icpu)
      call serxst (f(tsc)  , 'P3', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(f(tsc)  , 'P3', trnch,ni,heurser,1.,  -1, icpu)
      call serxst (f(tss)  , 'P4', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(f(tss)  , 'P4', trnch,ni,heurser,1.,  -1, icpu)
*
*     CONVECTIVE AND STRATIFORM PRECIPITATION RATES
      call serxst (v(rc)   , 'RC', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(v(rc)   , 'RC', trnch,ni,heurser,1.,  -1, icpu)
      call serxst (v(rr)   , 'RR', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(v(rr)   , 'RR', trnch,ni,heurser,1.,  -1, icpu)
*
*     ACCUMULATED PRECIPITATION
      call serxst (v(pr)   , 'PR', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(v(pr)   , 'PR', trnch,ni,heurser,1.,  -1, icpu)
      call serxst (v(pc)   , 'PC', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(v(pc)   , 'PC', trnch,ni,heurser,1.,  -1, icpu)
      call serxst (v(ae)   , 'AE', trnch,ni, 0.,    1.,  -1      )
      call mzonxst(v(ae)   , 'AE', trnch,ni,heurser,1.,  -1, icpu)
*
*     CLOUDS
      CALL SERXST (f(ccn)  , 'NU', trnch,ni, 0.,    1.,  -1      )
      CALL MZONXST(f(ccn)  , 'NU', trnch,ni,heurser,1.,  -1, ICPU)
      CALL SERXST (f(fn)   , 'NJ', trnch,ni, 0.,    1.,  -1      )
      CALL SERXST (f(qcbl) , 'QB', trnch,ni, 0.,    1.,  -1      )
*

      if ( iconvec.ge.3) then

*     NUAGES CONVECTIFS
      CALL SERXST (F(CCK)  , 'NC', trnch,ni, 0.,    1.,  -1      )
      CALL MZONXST(F(CCK)  , 'NC', trnch,ni,heurser,1.,  -1, ICPU)

      endif
*
*
*     EFFECT OF PRECIPITATION EVAPORATION
      do 28 k=1,nk-1
      do 28 i=1,ni
         work2d1(i,k)=min(0.,v(tcond+ik(i,k)))
28    continue
*
      call serxst (work2d1 , 'EP', trnch,ni,0.0,    1.       ,-1     )
      call mzonxst(work2d1 , 'EP', trnch,ni,heurser,d(pmoins),-2,icpu)
*
*
*
*     EXTRACTION OF DYNAMICS VARIABLES
*     --------------------------------
*
*     EXRACT TIME SERIES AND ZONAL DIAGNOSTICS ON NK LEVELS
      call mzoniv ( trnch, nk )
      call sersetm( 'KA', trnch, nk )
*
*
*     EXTRACTION OF WIND COMPONENTS
*
      CALL SERXST ( D(UPLUS) , 'UU' , TRNCH , NI , 0.0 , 1.0 , -1 )
      CALL SERXST ( D(VPLUS) , 'VV' , TRNCH , NI , 0.0 , 1.0 , -1 )
*
      CALL MVZNXST( D(UPLUS), D(VPLUS),'UU','VV',TRNCH,NI,1.0,-1,ICPU)
*
      CALL SERXST ( D(OMEGAP) , 'WW' , TRNCH , NI , 0.0 ,   1.0 , -1        )
      CALL MZONXST( D(OMEGAP) , 'WW' , TRNCH , NI , HEURE , 1.0 , -1 , ICPU )
*
*     TEMPERATURE
      CALL SERXST ( D(TPLUS), 'TT', TRNCH , NI , 0.0   , 1.0      , -1        )
      CALL MZONXST( D(TPLUS), 'TT', TRNCH , NI , HEURE , D(PPLUS) , -2 , ICPU )
*
*     POTENTIAL TEMPERATURE
      DO K=1,NK
*VDIR NODEP
         DO I=1,NI
            P   (I,K) = D(SIGM +IK(I,K))*D(PPLUS+IK(I,1))
            WORK2D1(I,K) = 1.E-5*P(I,K)
         END DO
      END DO
*
      CALL VSPOWN1 (WORK2D1,WORK2D1,-CAPPA,NI*NK)
*
      DO K=1,NK
         DO I=1,NI
            WORK2D1(I,K) = D(TPLUS+IK(I,K))*WORK2D1(I,K)
         END DO
      END DO
*
      CALL SERXST ( WORK2D1, 'TH', TRNCH , NI , 0.0   , 1.0 , -1        )
      CALL MZONXST( WORK2D1, 'TH', TRNCH , NI , HEURE , 1.0 , -1 , ICPU )
*
*     DEW POINT TEMPERATURE
*
*     FIND DEW POINT DEPRESSION FIRST (ES). SATURATION CALCULATED WITH 
*     RESPECT TO WATER ONLY (SINCE TD MAY BE COMPARED TO OBSERVED TEPHIGRAM).
      CALL MHUAES(WORK2D1,D(HUPLUS),D(TPLUS),-1.,P,3,.TRUE.,.FALSE.,NI,NK,NI)
*
      DO K=1,NK
*VDIR NODEP
         DO I=1,NI
            WORK2D1(I,K) = MIN(D(TPLUS+IK(I,K)),D(TPLUS+IK(I,K))-WORK2D1(I,K))
         END DO
      END DO
*
      CALL SERXST ( WORK2D1, 'TD', TRNCH , NI , 0.0   , 1.0 , -1        )
      CALL MZONXST( WORK2D1, 'TD', TRNCH , NI , HEURE , 1.0 , -1 , ICPU )
*
*     MOISTURE
      IF(WET) THEN
*
*        ELIMINATE NEGATIVE VALUES
*VDIR NODEP
         DO I=1,NI*NK
            WORK2D1(I,1) = MAX(0.0,D(HUPLUS+IK(I,1)))
         END DO
*
         CALL SERXST ( WORK2D1, 'HU', TRNCH, NI , 0.0,   1.0      , -1         )
         CALL MZONXST( WORK2D1, 'HU', TRNCH, NI , HEURE, D(PPLUS) , -2, ICPU  )
*
*        RELATIVE HUMIDITY
         MODP = 3
         tmp1 = 0.0
         IF (SATUCO) THEN
            DO K=1,NK
*VDIR NODEP
            DO I=1,NI
               WORK2D2(I,K) = D(SIGM+IK(I,K))*D(PPLUS+IK(I,1))

            END DO
            END DO
*
            CALL MFOHR(WORK2D1,D(HUPLUS), D(TPLUS),
     +                 TMP1, WORK2D2, MODP,NI,NK,NI)
*
         ELSE
            DO K=1,NK
            DO I=1,NI
               WORK2D2(I,K) = D(SIGM+IK(I,K))*D(PPLUS+IK(I,1))
            END DO
            END DO
*
            CALL MFOHRA(WORK2D1, D(HUPLUS), D(TPLUS),
     +                  TMP1, WORK2D2, MODP,NI,NK,NI)
*
         END IF
*
         CALL SERXST ( WORK2D1, 'HR', TRNCH, NI, 0.0   ,1.0, -1      )
         CALL MZONXST( WORK2D1, 'HR', TRNCH, NI, HEURE ,1.0, -1,ICPU)
*
         IF (ISTCOND. GE. 2) THEN
*
*           TOTAL CLOUD WATER
            CALL SERXST ( D(QCPLUS), 'QC', TRNCH, NI , 0.0,   1.0      , -1        )
            CALL MZONXST( D(QCPLUS), 'QC', TRNCH, NI , HEURE, D(PPLUS) , -2, ICPU )
*
            IF (ISTCOND. EQ. 5) THEN
*
*              SOLID CLOUD WATER AS SEEN BY RADIATION
               CALL SERXST ( F(IWC), 'QI', TRNCH, NI , 0.0,   1.0      , -1        )
               CALL MZONXST( F(IWC), 'QI', TRNCH, NI , HEURE, D(PPLUS) , -2, ICPU )
*
*              LIQUID CLOUD WATER AS SEEN BY RADIATION
               CALL SERXST ( F(LWC), 'QW', TRNCH, NI , 0.0,   1.0      , -1        )
               CALL MZONXST( F(LWC), 'QW', TRNCH, NI , HEURE, D(PPLUS) , -2, ICPU )
*
            ENDIF
*
         ENDIF
*
      ENDIF
*
*
*     SURFACE PRESSURE
      CALL SERXST ( D(PPLUS) , 'P0' , TRNCH , NI , 0.0  , 1.0 , -1         )
      CALL MZONXST( D(PPLUS) , 'P0' , TRNCH , NI, HEURE , 1.0 , -1 , ICPU )
*
*     MODULUS OF THE SURFACE WIND
      DO I = 1,NI
         WORK1D(I) = D(UPLUS+IK(I,NK))*D(UPLUS+IK(I,NK)) +
     %               D(VPLUS+IK(I,NK))*D(VPLUS+IK(I,NK))
      END DO
*
      call VSSQRT(WORK1D,WORK1D,NI)
*
      CALL SERXST ( WORK1D , 'VE' , TRNCH , NI , 0.0   , 1.0 , -1         )
      CALL MZONXST( WORK1D , 'VE' , TRNCH , NI , HEURE , 1.0 , -1 , ICPU )
*
*
      RETURN
      END