      SUBROUTINE OIFLUX11(QSENS,  TFLUX,  QEVAP,  EVAP,   QFLUX,        
     1                    CDH,    CDM,   SLIM,
     2                    ST,     SU,     SV,     SQ,     SRH,        
     3                    DRAG,   GC,     UGCM,   VGCM,   VMOD, 
     4                    TH,   
     5                    TA,     QA,     PRESSG, GT,     QG,              
     6                    ZDSLM,  ZDSLH,  ZREF,
     7                    ILG,    IL1,    IL2,    JL,     ISLFD        )   
C
C     * JUN 06/14 - M.Lazare.   New version for gcm18:
C     *                         - Add call to SLDIAG for ISLFD=1. This
C     *                           requires passing in ISLFD switch along
C     *                           with new input fields {ZDSLM,ZDSLH,ZREF}.
C     *                         - ZREF now passed, so no need to pass SGL,SHL,
C     *                           and is used in place of previous internal ZLEV.
C     *                         - Modification of screen-level calculations
C     *                           from Diana.
C     *                         - Replace data constants for
C     *                           {zero,one,two,ten} by actual values.
C     *                         - Correct calculation for SQ, ie same
C     *                           bugfix as in CLASST.
C     *                         - SICN removed and ice-only calculations
C     *                           for RATMOH (RATMOH=1.) with no SICN
C     *                           weighting. This is the proper approach
C     *                           for fractional water/ice.
C     * OCT 14/11 - M.LAZARE.   New version for gcm17:
C     *                         - Calls Diana's DRCOEF to be consistent
C     *                           with new class_v3.5, thus unused internal
C     *                           work arrays ZOMS,ZOHS removed.
C     * SEP 16/11 - M.LAZARE.   NEW VERSION FOR GCM15J:
C     *                         - Shortened CLASS4 common block consistent
C     *                           with new CLASS version.
C     * APR 23/10 - M.LAZARE.   PREVIOUS VERSION OIFLUX9 FOR GCM15I:
C     *                         - ADDS CALCULATION OF SCREEN RH THROUGH
C     *                           CALL TO NEW SUBROUTINE "SCREENRH".
C     *                           "SRH" IS THE DIAGNOSTIC OUTPUT PASSED
C     *                           BACK OUT.
C     * JAN 17/08 - M.LAZARE.   PREVIOUS VERSION OIFLUX8 FOR GCM15G/H:
C     *                         INCLUDE GUSTINESS EFFECTS FOR WINDSPEED
C     *                         THAT SHOULD HAVE BEEN INCORPORATED INTO
C     *                         THIS ROUTINE EARLIER ON (CONTAINED IN
C     *                         "VMODL" PASSED FROM CLASS DRIVER).
C     * JUN 15/06 - M.LAZARE.   PREVIOUS VERSION OIFLUX7 FOR GCM15F:
C     *                         - "ZERO","ONE","TWO","TEN","P9" DATA
C     *                            CONSTANTS ADDED AND USED IN CALLS
C     *                            TO INTRINSICS.
C     *                          - USES INTERNAL WORK ARRAYS INSTEAD
C     *                            OF PASSING IN "WRK" WORKSPACE FROM
C     *                            PHYSICS.  
C     *                          - REMOVE SETTING OF TRANSP TO ZERO -
C     *                            THIS IS DONE THROUGH INITIALIZATION
C     *                            ALREADY AT BEGINNING OF CLASS DRIVER.
C     *                            SINCE THIS IS THE ONLY PLACE WHERE
C     *                            TRANSP IS USED, IT IS REMOVED FROM
C     *                            THE ROUTINE ALLTOGETHER. 
C     * DEC 12/05 - K.VONSALZEN. PREVIOUS VERSION OIFLUX6 FOR GCM15D/E:
C     *                          - ADD GUSTINESS FIELD TO WIND SPEED 
C     *                            (PASSED IN FROM PHYSICS).
C     * DEC 14/04 - M.LAZARE. PREVIOUS VERSION OIFLUX5 FOR GCM15C:
C     *                       CALLS NEW DRCOEF5 VERSION.
C     * FEB 18/03 - M.LAZARE. PREVIOUS VERSION OIFLUX4 FOR GCM15B.
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C                                                                            
C     * INPUT/OUTPUT ARRAYS.                                            
C                                                                             
      REAL  , DIMENSION(ILG)  :: QSENS,TFLUX,QEVAP,EVAP,QFLUX,          
     1                           CDH,CDM,SLIM,ST,SU,SV,SQ,SRH,DRAG
C                                                                             
C     * INPUT ARRAYS.                                                         
C                                                                             
      REAL  , DIMENSION(ILG)  :: GC,UGCM,VGCM,VMOD,TH,
     1                           TA,QA,PRESSG,GT,QG,
     2                           ZDSLM,ZDSLH,ZREF
C                                                                             
C     * INTERNAL WORK ARRAYS.
C                                                                             
      REAL  , DIMENSION(ILG)  :: RIB,GTV,TAV,CDO,CLHTOP,RATMOH,
     1                           VA,ZOSCLM,ZOSCLH,
     2                           CRIB,CFLUX,FCT,ZRSLDM,ZRSLDH,
     3                           FRAC,ZOM,ZOH
C
      INTEGER ITER   (ILG)
C
      COMMON /CLASS1/ DELT,TFREZ                                               
      COMMON /CLASS2/ RGAS,RGASV,GRAV,SBC,VKC,CT,VMIN
      COMMON /CLASS4/ HCPW,HCPICE,HCPSOL,HCPOM,HCPSND,HCPCLY,
     1                SPHW,SPHICE,SPHVEG,SPHAIR,RHOW,RHOICE,
     2                TCGLAC,CLHMLT,CLHVAP
C---------------------------------------------------------------------
C     * DETERMINE NEUTRAL DRAG PARAMETERS.
C     * ASSUME NEUTRAL DRAG COEFFICIENT FOR HEAT OVER SEA-ICE IS
C     * 2.5X THAT OF OPEN WATER.
C                                                                              
      DO 100 I=IL1,IL2                                                        
          VA(I)=MAX(VMIN,VMOD(I))
          CRIB(I)=-GRAV*ZREF(I)/(TH(I)*(1.+0.61*QA(I))*(VA(I)**2))
          IF(GC(I).GT.-0.5.AND.GC(I).LT.0.5) THEN                   
              FM2=(1.+SQRT(CT)*LOG(ZREF(I)/10.)/VKC)**2             
              CDO(I)=CT/FM2                                         
C             QUAD=(MAX(VA(I)-10. , 0.))**2                         
              QUAD=0.                                               
              RATMOH(I)=0.001*MAX(0.9,0.61+0.063*VA(I)+0.05*QUAD)/  
     1                  (CDO(I)*FM2)                                
              CLHTOP(I)=CLHVAP                                      
          ELSE IF(GC(I).GT.0.5) THEN                                   
              CTEFF=2.5*CT
              FM2=(1.+SQRT(CTEFF)*LOG(ZREF(I)/10.)/VKC)**2 
              CDO(I)=CTEFF/FM2                          
              RATMOH(I)=1.
              CLHTOP(I)=CLHVAP+CLHMLT                               
          ENDIF                                                     
          DRAG(I)=CDO(I)
          GTV(I)=GT(I)*(1.0+0.61*QG(I))                                      
          TAV(I)=TA(I)*(1.0+0.61*QA(I))
          ZOSCLH(I)=EXP(-VKC*SQRT(CDO(I)*RATMOH(I))/CDO(I)) 
          ZOSCLM(I)=EXP(-VKC/SQRT(CDO(I)*RATMOH(I)))
          FCT(I)=1.
          ITER(I)=1
  100 CONTINUE                                                               
C
C     * CALCULATE SURFACE DRAG COEFFICIENTS (STABILITY-DEPENDANT) AND
C     * OTHER RELATED QUANTITIES.
C     * NOTE IN THE ABOVE LOOP THAT FCT=1. (WHOLE GRID CELL OPEN
C     * WATER) AND ITER=1 (ALWAYS DOING LOOP, ONE TIME ONLY).
C
      CALL DRCOEF(CDM,CDH,RIB,CFLUX,QG,QA,ZOSCLM,ZOSCLH,
     1            CRIB,GTV,TAV,VA,FCT,ITER,
     2            ILG,IL1,IL2)
C
C     * REMAINING CALCULATIONS.
C
      DO 500 I=IL1,IL2
          QLWOUT=SBC*GT(I)*GT(I)*GT(I)*GT(I)                                  
          EA=QA(I)*PRESSG(I)/(0.622+0.378*QA(I))                              
          PADRY=PRESSG(I)-EA                                                  
          RHOAIR=PADRY/(RGAS*TA(I))+EA/(RGASV*TA(I))                          
          SLIM(I)=RHOAIR*CFLUX(I)
          if (abs(SLIM(I)) > 1e10) then
            write(6,*)"oiflux11:  SLIM=",SLIM(I),
     1        "  CFLUX=",CFLUX(I),"  rhoair=",rhoair
          endif
          QSENS(I)=RHOAIR*SPHAIR*CFLUX(I)*(GT(I)-TA(I))                       
          IF(TA(I).GE.TFREZ) THEN                                              
              A=17.269                                                         
              B=35.86                                                          
          ELSE                                                                 
              A=21.874                                                         
              B=7.66                                                           
          ENDIF                                                                
          WA=QA(I)/(1.0-QA(I))                                                 
          CONST=LOG(WA*PADRY/(0.622*611.0))                                   
          TADP=(B*CONST-A*TFREZ)/(CONST-A)                                     
          EVAP(I)=RHOAIR*CFLUX(I)*(QG(I)-QA(I))                               
          IF(EVAP(I).LT.0.0.AND.GT(I).GE.TADP) EVAP(I)=0.0                    
          QEVAP(I)=CLHTOP(I)*EVAP(I)                                           
          TFLUX(I)=-QSENS(I)/(RHOAIR*SPHAIR)                                  
          QFLUX(I)=-EVAP(I)/RHOAIR                                            
          FRAC (I)=1.
          ZOM(I)=ZOSCLM(I)*ZREF(I)
          ZOH(I)=ZOSCLH(I)*ZREF(I)
          ZRSLDM(I)=ZREF(I)+ZOM(I)
          ZRSLDH(I)=ZREF(I)+ZOM(I)
  500 CONTINUE
C
C     * SCREEN LEVEL CALCULATIONS.
C
      IF(ISLFD.EQ.0)                                         THEN
        DO 550 I=IL1,IL2
          FACTM=ZDSLM(I)+ZOM(I)
          FACTH=ZDSLH(I)+ZOM(I)
          RATIOM=SQRT(CDM(I))*LOG(FACTM/ZOM(I))/VKC
          RATIOM=MIN(RATIOM,1.)
          RATIOH=SQRT(CDM(I))*LOG(FACTH/ZOH(I))/VKC 
          RATIOH=MIN(RATIOH,1.)
          IF(RIB(I).LT.0.)  THEN
             RATIOH=RATIOH*CDH(I)/CDM(I)
             RATIOH=MIN(RATIOH,(FACTH/ZRSLDH(I))**(1./3.)) 
          ENDIF                                
          ST(I)=GT(I)-(MIN(RATIOH,1.))*(GT(I)-TA(I))
          SQ(I)=QG(I)-(MIN(RATIOH,1.))*(QG(I)-QA(I))
          SU(I)=RATIOM*UGCM(I)
          SV(I)=RATIOM*VGCM(I)
  550   CONTINUE                                                              
      ELSE IF(ISLFD.EQ.1)                                    THEN
        CALL SLDIAG(SU,SV,ST,SQ,
     1              CDM,CDH,UGCM,VGCM,TA,QA,
     2              GT,QG,ZOM,ZOH,FRAC,ZRSLDM,
     3              ZDSLM,ZDSLH,ILG,IL1,IL2,JL)
      ELSE
        CALL XIT('OIFLUX11',-1)
      ENDIF
C
C     * CALCULATE SCREEN RH.
C
      CALL SCREENRH(SRH,ST,SQ,PRESSG,FRAC,ILG,IL1,IL2)
C
      RETURN                                                                
      END                                                                      
