      PROGRAM DPDETA                                                                                                                
C     PROGRAM DPDETA(XIN,      LNSP,      XOUT,      INPUT,      OUTPUT,)       J2                                                  
C    1         TAPE1=XIN,TAPE2=LNSP,TAPE3=XOUT,TAPE5=INPUT,TAPE6=OUTPUT)                                                            
C     ------------------------------------------------------------------        J2                                                  
C                                                                               J2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       J2                                                  
C     DEC 06/93 - D.LIU                                                                                                             
C                                                                               J2                                                  
CDPDETA  - DERIVATIVE OF PRESSURE WITH RESPECT TO ETA COORDINATES       2  1 C  J1                                                  
C                                                                               J2                                                  
CAUTHOR  - D. LIU                                                               J3                                                  
C                                                                               J3                                                  
CPURPOSE - COMPUTES DERIVATIVE OF PRESSURE WITH RESPECT TO ETA COORDINATES      J3                                                  
C                                                                               J3                                                  
CINPUT FILES...                                                                 J3                                                  
C                                                                               J3                                                  
C      XIN  = INPUT SERIES OF ETA LEVEL GRID FIELDS .                           J3                                                  
C      LNSP = INPUT SERIES OF LN(SURFACE PRESSURE) IN MB.                       J3                                                  
C                                                                               J3                                                  
COUTPUT FILE...                                                                 J3                                                  
C                                                                               J3                                                  
C      XOUT = DP/D(ETA)  (MB)                                                   J3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               J5                                                  
C      LEVTYP = FULL FOR MOMENTUM VARIABLE, AND                                 J5                                                  
C               HALF FOR THERMODYNAMIC ONE.                                     J5                                                  
C      LAY    = DEFINES THE POSITION OF LAYER INTERFACES IN RELATION            J5                                                  
C               TO LAYER CENTRES (SEE BASCAL).                                  J5                                                  
C               (ZERO DEFAULTS TO THE FORMER STAGGERING CONVENTION).            J5                                                  
C      ICOORD = 4H ETA/4HET15 FOR ETA LEVELS OF INPUT FIELDS.                   J5                                                  
C      SIGTOP = VALUE OF SIGMA AT TOP OF DOMAIN FOR VERTICAL INTEGRAL.          J5                                                  
C               IF .LT.0., THEN INTERNALLY DEFINED BASED ON LEVTYP              J5                                                  
C               FOR UPWARD COMPATIBILITY.                                       J5                                                  
C      PTOIT  = PRESSURE (PA) OF THE RIGID LID OF THE MODEL.                    J5                                                  
C                                                                               J5                                                  
C      NOTE - LAY AND LEVTYP DEFINE THE TYPE OF LEVELLING FOR THE VARIABLE.     J5                                                  
C                                                                               J5                                                  
CEXAMPLE OF INPUT CARD...                                                       J5                                                  
C                                                                               J5                                                  
C*DPDETA   HALF              3 ET15-1.00       50.                              J5                                                  
C------------------------------------------------------------------------------                                                     
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK,SPEC                                                                                                               
                                                                                                                                    
      INTEGER LEV(100)                                                                                                              
                                                                                                                                    
      REAL ETA (100),ETAB (100),E (2)                                                                                               
      REAL AB  (100),BB   (100)                                                                                                     
                                                                                                                                    
      COMMON /LEVELS/ SIGB  (18528,2)                                                                                               
      COMMON /BLANCK/ PSMBLN(18528), GG(18528)                                                                                      
      REAL             PRESS(18528)                                                                                                 
      EQUIVALENCE     (PRESS,PSMBLN)                                                                                                
                                                                                                                                    
      COMMON /ICOM  /IBUF(8),IDAT(18528)                                                                                            
                                                                                                                                    
      DATA MAXX/18528/, MAXLEV/100/                                                                                                 
                                                                                                                                    
      PARAMETER (P0=101320.E0)                                                                                                      
C---------------------------------------------------------------------                                                              
      NFIL = 5                                                                                                                      
      CALL JCLPNT(NFIL,1,2,3,5,6)                                                                                                   
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      REWIND 3                                                                                                                      
                                                                                                                                    
C     * READ CONTROL DIRECTIVES.                                                                                                    
                                                                                                                                    
      READ(5,5010,END=900) LEVTYP,LAY,ICOORD,SIGTOP,PTOIT                       J4                                                  
      IF(ICOORD.EQ.NC4TO8("    ")) ICOORD=NC4TO8(" SIG")                                                                            
      IF(ICOORD.EQ.NC4TO8(" SIG")) THEN                                                                                             
        PTOIT=MAX(PTOIT,0.00E0)                                                                                                     
      ELSE                                                                                                                          
        PTOIT=MAX(PTOIT, 0.100E-09)                                                                                                 
      ENDIF                                                                                                                         
      WRITE(6,6005) LEVTYP,LAY,ICOORD,SIGTOP,PTOIT                                                                                  
                                                                                                                                    
C     * GET ETA VALUES OF LAYER CENTRES FROM XIN.                                                                                   
                                                                                                                                    
      CALL FILEV (LEV,ILEV,IBUF,1)                                                                                                  
      IF(ILEV.LT.1 .OR. ILEV.GT.MAXLEV)  CALL      XIT('DPDETA',-1)                                                                 
      WRITE(6,6007) ILEV,(LEV(L),L=1,ILEV)                                                                                          
                                                                                                                                    
      CALL LVDCODE(ETA,LEV,ILEV)                                                                                                    
      DO 150 L=1,ILEV                                                                                                               
         ETA(L)=ETA(L)*0.001E0                                                                                                      
  150 CONTINUE                                                                                                                      
                                                                                                                                    
C     * DETERMINE THE FIELD SIZE.                                                                                                   
                                                                                                                                    
      KIND=IBUF(1)                                                                                                                  
      NAME=IBUF(3)                                                                                                                  
      SPEC=(KIND.EQ.NC4TO8("SPEC") .OR. KIND.EQ.NC4TO8("FOUR"))                                                                     
      NWDS=IBUF(5)*IBUF(6)                                                                                                          
      IF(SPEC) CALL                                XIT('DPDETA',-2)                                                                 
                                                                                                                                    
                                                                                                                                    
C     * EVALUATE LAYER INTERFACES FROM LEVTYP AND LAY.                                                                              
                                                                                                                                    
      IF    (LEVTYP.EQ.NC4TO8("FULL")) THEN                                                                                         
         CALL BASCAL (ETAB,IBUF, ETA,ETA,ILEV,LAY)                                                                                  
      ELSEIF(LEVTYP.EQ.NC4TO8("HALF")) THEN                                                                                         
         CALL BASCAL (IBUF,ETAB, ETA,ETA,ILEV,LAY)                                                                                  
      ELSE                                                                                                                          
         CALL                                      XIT('DPDETA',-3)                                                                 
      ENDIF                                                                                                                         
                                                                                                                                    
                                                                                                                                    
C     * EVALUATE THE PARAMETERS A AND B OF THE VERTICAL DISCRETIZATION                                                              
C     * FOR THE LAYER BASES.                                                                                                        
                                                                                                                                    
      CALL COORDAB (AB,BB, ILEV,ETAB,ICOORD,PTOIT)                                                                                  
                                                                                                                                    
C---------------------------------------------------------------------                                                              
      NSETS=0                                                                                                                       
                                                                                                                                    
  200 CONTINUE                                                                                                                      
                                                                                                                                    
      CALL GETFLD2 (2,PSMBLN,KIND,-1,NC4TO8("LNSP"),1,IBUF,MAXX,OK)                                                                 
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      IF(.NOT.OK) THEN                                                                                                              
         IF(NSETS.GT.0) THEN                                                                                                        
            WRITE(6,6025) IBUF                                                                                                      
            WRITE(6,6020) NSETS                                                                                                     
            CALL                                   XIT('DPDETA',0)                                                                  
         ELSE                                                                                                                       
            CALL                                   XIT('DPDETA',-4)                                                                 
         ENDIF                                                                                                                      
      ENDIF                                                                                                                         
      ITIM=IBUF(2)                                                                                                                  
      DO 300 I=1,NWDS                                                                                                               
         PRESS(I)=100.0E0*EXP(PSMBLN(I))                                                                                            
  300    CONTINUE                                                                                                                   
                                                                                                                                    
C     * DEFINE TOP SIGMA VALUE.                                                                                                     
                                                                                                                                    
      IGH = 1                                                                                                                       
      LOW = 2                                                                                                                       
      IF(SIGTOP.LT.0.E0)THEN                                                                                                        
         IF(LEVTYP.EQ.NC4TO8("FULL"))THEN                                                                                           
            DO 440 I=1,NWDS                                                                                                         
               SIGB(I,IGH)=0.0E0                                                                                                    
  440       CONTINUE                                                                                                                
         ELSE                                                                                                                       
            DO 450 I=1,NWDS                                                                                                         
               SIGB(I,IGH)=PTOIT/PRESS(I)                                                                                           
  450       CONTINUE                                                                                                                
         ENDIF                                                                                                                      
      ELSE                                                                                                                          
         DO 460 I=1,NWDS                                                                                                            
            SIGB(I,IGH)=SIGTOP                                                                                                      
  460    CONTINUE                                                                                                                   
      ENDIF                                                                                                                         
                                                                                                                                    
      E(IGH)=PTOIT/P0                                                                                                               
                                                                                                                                    
C     * PROCESS ALL LEVELS. SIGMA MULTIPLIED BY PRESS TO GIVE PRESSURE (MB).                                                        
                                                                                                                                    
      DO 800 L=1,ILEV                                                                                                               
                                                                                                                                    
         CALL GETFLD2 (1,GG,KIND,ITIM,NAME,LEV(L),IBUF,MAXX,OK)                                                                     
         IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                          
         IF(.NOT.OK)  CALL                         XIT('DPDETA',-10-L)                                                              
                                                                                                                                    
         CALL NIVCAL(SIGB(1,LOW), AB(L),BB(L),PRESS,1,NWDS,NWDS)                                                                    
         E(LOW)=ETAB(L)                                                                                                             
                                                                                                                                    
         DO 500 I=1,NWDS                                                                                                            
            GG(I) = (SIGB(I,LOW)-SIGB(I,IGH))*PRESS(I)/100.E0                                                                       
     1              /(E(LOW)-E(IGH))                                                                                                
  500    CONTINUE                                                                                                                   
                                                                                                                                    
         KEEP= LOW                                                                                                                  
         LOW = IGH                                                                                                                  
         IGH = KEEP                                                                                                                 
                                                                                                                                    
                                                                                                                                    
C     * PUT THE RESULT ONTO FILE XOUT.                                                                                              
                                                                                                                                    
      IBUF(3)=NC4TO8("DPDE")                                                                                                        
      CALL PUTFLD2 (3,GG,IBUF,MAXX)                                                                                                 
      NSETS=NSETS+1                                                                                                                 
                                                                                                                                    
  800 CONTINUE                                                                                                                      
                                                                                                                                    
      GO TO 200                                                                                                                     
                                                                                                                                    
C     * E.O.F. ON INPUT.                                                                                                            
                                                                                                                                    
  900 CALL                                         XIT('DPDETA',-5)                                                                 
C---------------------------------------------------------------------                                                              
 5010 FORMAT(11X,A4,I15,1X,A4,E5.0,E10.0)                                       J4                                                  
 6005 FORMAT('  LEVTYP = ',A4,                                                                                                      
     1       ', LAY = ',0P,I5,', COORD=',1X,A4,', SIGTOP =',F15.10,                                                                 
     2       ', P.LID (PA)=',E10.3)                                                                                                 
 6007 FORMAT(' ILEV,ETA =',I5,2X,20I5/(18X,20I5))                                                                                   
 6020 FORMAT('0',I6,' SETS WERE PROCESSED')                                                                                         
 6025 FORMAT(' ',A4,I10,2X,A4,5I10)                                                                                                 
      END                                                                                                                           
