C                                                                       
      SUBROUTINE MDIV(J1,JHR,FD,FD1,FD2,FD3,JJ,II,X,ICND,*)
C     FEB 1978    GLAHN       MDL    FOR M201                           
C     JAN 1985    ERICKSON    MDL    (REVISED)                          
C     JAN 1991    DALLAVALLE   MDL   (REVISED TO USE REL. HUM., INSTEAD 
C                 OF DEW POINT, FOR COMPUTATION OF THE MOISTURE DIV.    
C                 FROM THE NGM.  NOTE THAT REL. HUM. IS THE MOISTURE    
C                 VARIABLE IN THE NGM.  CHANGES WERE ALSO MADE TO       
C                 ALLOW THE COMPUTATION AT 950, 850, 700, 500, OR 300 MB
C                 IN THE NGM.)                                          
C     JUN 1991    GILBERT      MDL    FORTRAN 77                        
C     PURPOSE                                                           
C        TO COMPUTE MOISTURE DIVERGENCE IN BOUNDARY LAYER OR AT 850,    
C        700, OR 500 MB.   ONLY THESE LEVELS ARE AVAILABLE IN THE       
C        LFM ARCHIVE.                                                   
C     VARIABLES                                                         
C              J1 = FIRST WORD OF MDL ID OF FORM YYYXXX.                
C                   IF XXX=114HEX, CORRECT COMPUTATION IS MADE USING    
C                   GRID-ORIENTED WINDS. THE INCORRECT COMPUTATION      
C                   USING EARTH-ORIENTED WINDS( XXX=113HEX) IS NO LONGER
C                   ALLOWED. (1985)                                     
C             JHR = TAU CORRESPONDING TO J1.                            
C       FD(JJ,II) = MOISTURE DIVERGENCE RETURNED.                       
C      FD1(JJ,II) = WORK ARRAY.                                         
C      FD2(JJ,II) = WORK ARRAY.                                         
C      FD3(JJ,II) = WORK ARRAY.                                         
C            ICND = 0 FOR GOOD RETURN.                                  
C               * = RETURN TO STATEMENT NUMBER.                         
      DIMENSION FD(JJ,II),FD1(JJ,II),FD2(JJ,II),FD3(JJ,II),X(1)
      COMMON/MAPSCL/POL2EQ(8),XMESHM(8)
      COMMON/M90201/JDATE,NCYCLE,MODP1,IG,JUNPK(8,4,7)
      DATA MASK/X'FFF'/
      DATA ITHETA/X'225023'/,IRH/X'225130'/,ISFC/X'231010'/
      DATA I060/X'060'/,I225/X'225'/,I160/X'160'/,I205/X'205'/,
     1     I021/X'021'/,I114/X'114'/,I226/X'226'/,I130/X'130'/,
     1     I020/X'020'/,I120/X'120'/
      K1=ISHFT(J1,-12)
C        K1=LEVEL                                                       
      K2=IAND(J1,MASK)
C        K2=TYPE (SHOULD BE 114)                                        
      GRD=POL2EQ(MODP1)
      DG=XMESHM(MODP1)
      IX=JUNPK(MODP1,IG,6)
      IY=JUNPK(MODP1,IG,7)
      DGX2=2.*DG
      III=II-1
      JJJ=JJ-1
      DO 120 K=1,II
      DXSQ=(IX-K)**2
      DO 120 J=1,JJ
      DYSQ=(IY-J)**2
      FD3(J,K)=(GRD-(DXSQ+DYSQ))/(GRD+(DXSQ+DYSQ))
      FD3(J,K)=1.866/(1.+FD3(J,K))
 120  CONTINUE
C        MAP SCALE FACTORS IN FD3( , ).                                 
C        CHECK IF THIS IS THE NGM (MODEL ID PLUS 1 EQUALS 7)            
      IF(MODP1.EQ.7)GO TO 1400
C        OTHERWISE, IT MUST BE THE LFM OR GLOBAL MODEL                  
      IF(K1.EQ.I225)GO TO 128
C        MOISTURE DIVERGENCE AT A LEVEL OTHER THAN BOUNDARY LAYER.      
C        USE DEW POINT RATHER THAN RELATIVE HUMIDITY.                   
      IDP=IOR(ISHFT(K1,12),I021)
C        IDP=IDENTIFIER FOR DEW POINT AT LEVEL INDICATED IN J1.         
      CALL GTPRED(IDP,JHR,FD1,JJ,II,X,3,X,ICND,*240)
C        DEW POINT IN FD1( , ).                                         
      PRES=0.
      IF(K1.EQ.I205)PRES=850.
      IF(K1.EQ.I160)PRES=700.
      IF(K1.EQ.I120)PRES=500.
      IF(PRES.NE.0.)GO TO 122
C        NOT ONE OF LEVELS 850, 700, OR 500 MB.                         
      ICND=1
      GO TO 240
 122  DO 123 K=1,II
      DO 123 J=1,JJ
      TEMP=FD1(J,K)-273.16
      RATIO=(17.269*TEMP)/(TEMP+237.3)
      ES=6.1078*EXP(RATIO)
      FD(J,K)=(.622*ES)/(PRES-.378*ES)
 123  CONTINUE
      GO TO 141
 128  CALL GTPRED(ITHETA,JHR,FD1,JJ,II,X,3,X,ICND,*240)
C        B.L. POTENTIAL TEMP IN FD1( , ).                               
      CALL GTPRED(IRH,JHR,FD2,JJ,II,X,3,X,ICND,*240)
C        B.L.  REL. HUM. IN FD2( , ).                                   
      CALL GTPRED(ISFC,JHR,FD,JJ,II,X,3,X,ICND,*240)
C     SFC PRESSURE IN FD( , ).                                          
      DO 130 K=1,II
      DO 130 J=1,JJ
      FD1(J,K)=FD1(J,K)*((FD(J,K)-25.)/1000.)**.286
 130  CONTINUE
C        CALCULATE SATURATED VAPOR PRESSURE AND SATURATED SPECIFIC      
C        HUMIDITY FROM TEMPERATURE  (STACKPOLE'S METHOD).               
      DO 140 K=1,II
      DO 140 J=1,JJ
      TEMP=FD1(J,K)-273.16
      RATIO=(17.269*TEMP)/(TEMP+237.3)
      ES=6.1078*EXP(RATIO)
      PRES=FD(J,K)-25.
      FD1(J,K)=(.622*ES)/(PRES-.378*ES)
      FD(J,K)=FD1(J,K)*FD2(J,K)*.01
 140  CONTINUE
C     SPECIFIC HUMIDITY IN FD( , ).                                     
      GO TO 141
 1400 PRES=0.
      IF(K1.EQ.I226)PRES=950.
      IF(K1.EQ.I205)PRES=850.
      IF(K1.EQ.I160)PRES=700.
      IF(K1.EQ.I120)PRES=500.
      IF(K1.EQ.I060)PRES=300.
C        IF THE PRESSURE STILL EQUALS 0, THEN THE MOISTURE DIVERGENCE   
C        CAN NOT BE COMPUTED AT THE LEVEL INDICATED BY THE ID.          
      IF(PRES.NE.0)GO TO 1410
      ICND=1
      GO TO 240
C        IDENTIFY THE RELATIVE HUMIDITY AT THE LEVEL DESIRED            
 1410 IRELH=IOR(ISHFT(K1,12),I130)
C        NOW OBTAIN THE RELATIVE HUMIDITY AND STORE THE VALUE IN FD2    
      CALL GTPRED(IRELH,JHR,FD2,JJ,II,X,3,X,ICND,*240)
C        IDENTIFY THE TEMPERATURE AT THE LEVEL DESIRED                  
      ITMP=IOR(ISHFT(K1,12),I020)
C        NOW OBTAIN THE TEMPERATURE AND STORE THE VALUE IN FD1          
      CALL GTPRED(ITMP,JHR,FD1,JJ,II,X,3,X,ICND,*240)
C        CALCULATE SATURATED VAPOR PRESSURE AND SATURATED SPECIFIC      
C        HUMIDITY FROM TEMPERATURE  (STACKPOLE'S METHOD).               
      DO 1420 K=1,II
      DO 1420 J=1,JJ
      TEMP=FD1(J,K)-273.16
      RATIO=(17.269*TEMP)/(TEMP+237.3)
      ES=6.1078*EXP(RATIO)
      FD1(J,K)=(.622*ES)/(PRES-.378*ES)
C        COMPUTE THE SPECIFIC HUMIDITY FROM THE RELATIVE HUMIDITY       
      FD(J,K)=FD1(J,K)*FD2(J,K)*.01
 1420 CONTINUE
C        THE SPECIFIC HUMIDITY IS LOCATED IN ARRAY FD                   
 141  IWND=IOR(ISHFT(K1,12),I060)
      CALL GTPRED(IWND,JHR,FD1,JJ,II,X,3,X,ICND,*240)
C        B.L.  U-WIND IN FD1( , ).                                      
      CALL GTPRED(IWND+1,JHR,FD2,JJ,II,X,3,X,ICND,*240)
C        B.L.  V-WIND IN FD2( , ).                                      
      IF(K2.EQ.I114)GO TO 155
C        INCORRECT J1                                                   
      ICND=1
      GO TO 240
 155  DO 160 K=1,II
      DO 160 J=1,JJ
      FD1(J,K)=FD1(J,K)*FD(J,K)
      FD2(J,K)=FD2(J,K)*FD(J,K)
 160  CONTINUE
C        PRODUCT OF WIND COMPONENTS AND SPECIFIC HUMIDITY IN            
C        FD1( , ) AND FD2( , ).                                         
      DO 170 K=2,III
      DO 170 J=2,JJJ
      FD(J,K)=(FD3(J,K)/DGX2)*(FD1(J,K+1)-FD1(J,K-1)+FD2(J+1,K)-
     1 FD2(J-1,K))
 170  CONTINUE
      DO 180 J=1,JJ
      FD(J,1)=(FD3(J,1)/DG)*(FD1(J,2)-FD1(J,1))
      FD(J,II)=(FD3(J,II)/DG)*(FD1(J,II)-FD1(J,III))
 180  CONTINUE
      DO 190 K=2,III
      FD(1,K)=(FD3(1,K)/DGX2)*(FD1(1,K+1)-FD1(1,K-1))
      FD(JJ,K)=(FD3(JJ,K)/DGX2)*(FD1(JJ,K+1)-FD1(JJ,K-1))
 190  CONTINUE
      DO 200 K=1,II
      FD(1,K)=FD(1,K)+(FD3(1,K)/DG)*(FD2(2,K)-FD2(1,K))
      FD(JJ,K)=FD(JJ,K)+(FD3(JJ,K)/DG)*(FD2(JJ,K)-FD2(JJJ,K))
 200  CONTINUE
      DO 210 J=2,JJJ
      FD(J,1)=FD(J,1)+(FD3(J,1)/DGX2)*(FD2(J+1,1)-FD2(J-1,1))
      FD(J,II)=FD(J,II)+(FD3(J,II)/DGX2)*(FD2(J+1,II)-FD2(J-1,II))
 210  CONTINUE
      DO 230 I=1,II
      DO 230 J=1,JJ
      FD(J,I)=FD(J,I)*1.0E+8
 230  CONTINUE
      ICND=0
 240  RETURN 1
      END
