!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
! -----------------------------------------------
! MSC-RPN COMM Group Licence/DISCLAIMER version 1
! Last Modified: Septembre 1st, 2001
! 
! MC2 : Canadian Mesoscale Compressible Model
!       Division de la Recherche en Prevision Numerique
!       Environnement Canada
! 
! This is a free software;
! you can redistribute it and/or modify it under the terms
! of the MSC-RPN COMM Group Licence/Disclaimer that can be found
! at http://www.cmc.ec.gc.ca/rpn/modcom/licence.html ; either
! version 1 or (at your option) any later version.
! 
! This software 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 MSC-RPN COMM Group Licence/Disclaimer for more details.
! 
! You should have received a copy of the MSC-RPN COMM Group
! Licence/Disclaimer along with this software; if not, write to the
! MSC-RPN COMM Group, 2121 TransCanada, suite 564, Dorval (Quebec),
! CANADA, H9P 1J3 ; or send a mail to comm.ec.gc.ca
! -----------------------------------------------

*** S/P WATER
*
#include "phy_macros_f.h"
      SUBROUTINE WATER ( bus, bussiz,
     $                   ptsurf, ptsurfsiz,
     $                   trnch, kount,
     $                   n, m, nk,seaspray,charnockl)
*
#include "impnone.cdk"
#include "partopo.cdk"
*
      integer bussiz, kount, trnch
      real bus(bussiz)
      integer ptsurfsiz
      integer ptsurf(ptsurfsiz)
*
*
      INTEGER N, M, NK
*
*Author
*          J. Mailhot, S. Belair and B. Bilodeau (Dec 1998)
*
*Revisions
* 001      B. Bilodeau (Nov 2000) - New comdeck sfcbus.cdk
* 002      B. Bilodeau (Jan 2001) - Automatic arrays
* 003      B. Bilodeau (Aug 2001) - LOCBUS
* 004      B. Bilodeau (Feb 2002) - Add Z0TCST option
*
*Object
*          Calculate: - surface roughness length (Z0) over open water
*                       (not covered by ice) using Charnock's relation with a 
*                       BETA parameter (optional sea state dependency);
*                     - surface fluxes of momentum, heat and moisture.
*
*Arguments
*
*             - Input/Output -
* BUS         Bus for the WATER surface scheme
*
*             - Input -
* BUSSIZ      dimension of bus
* PTSURF      surface pointers
* PTSURFSIZ   dimension of ptsurf
* TRNCH       row number
* KOUNT       timestep number
* N           horizontal dimension (row length)
* M           horizontal dimensions of fields
*             (not used for the moment)
* NK          vertical dimension
*
*
*Notes
*          Z0 = BETA*USTAR**2/GRAV (in metres) with minimum value
*          Z0MIN and a maximum value Z0MAX
*
*IMPLICITES
*
#include "consphy.cdk"
*
#include "clefcon.cdk"
#include "options.cdk"
#include "couplecp.cdk"
*
      integer ptr, x
      INTEGER I, J, K ,IC, JC,IW,JW
*
*
*
*MODULES
*
      EXTERNAL FLXSURF3,FLXSURF2
*
*
****************************************************
*     AUTOMATIC ARRAYS
****************************************************
*
      AUTOMATIC ( VMOD  , REAL , (N) )
      AUTOMATIC ( SCR1  , REAL , (N) )
      AUTOMATIC ( SCR2  , REAL , (N) )
      AUTOMATIC ( SCR3  , REAL , (N) )
      AUTOMATIC ( SCR4  , REAL , (N) )
      AUTOMATIC ( SCR5  , REAL , (N) )
      AUTOMATIC ( SCR6  , REAL , (N) )
      AUTOMATIC ( SCR7  , REAL , (N) )
      AUTOMATIC ( ZTN   , REAL , (N) )
      AUTOMATIC ( ZUN   , REAL , (N) )
      AUTOMATIC ( CP    , REAL , (N) )
      AUTOMATIC ( CHW3  , REAL , (N) )
*
****************************************************
*
      REAL BETA, RHO
*
      REAL ALVIS_WAT, CMU, CTU, FC_WAT
      REAL FV_WAT
      REAL HST_WAT, HU, ILMO_WAT
      REAL PS, QS, TH, TS, TT, UU, VV
      REAL Z0H, Z0M
      REAL ZALFAQ, ZALFAT, ZDLAT, ZFCOR
      REAL ZFTEMP, ZFVAP, ZQDIAG, ZTDIAG
      REAL ZTSURF, ZTSRAD, ZUDIAG, ZVDIAG, ZUE2, ZZA
CZWQ
      REAL ZWW3LAT,ZWW3LON
*
      POINTER (IALVIS_WAT , ALVIS_WAT  (1) )
      POINTER (ICMU       , CMU        (1) )
      POINTER (ICTU       , CTU        (1) )
      POINTER (IFC   _WAT , FC   _WAT  (1) )
      POINTER (IFV   _WAT , FV   _WAT  (1) )
      POINTER (IHST  _WAT , HST  _WAT  (1) )
      POINTER (IHU        , HU         (1) )
      POINTER (IILMO _WAT , ILMO _WAT  (1) )
      POINTER (IPS        , PS         (1) )
      POINTER (IQS        , QS         (1) )
      POINTER (ITH        , TH         (1) )
      POINTER (ITS        , TS         (1) )
      POINTER (ITT        , TT         (1) )
      POINTER (IUU        , UU         (1) )
      POINTER (IVV        , VV         (1) )
      POINTER (IZ0H       , Z0H        (1) )
      POINTER (IZ0M       , Z0M        (1) )
      POINTER (IZALFAQ    , ZALFAQ     (1) )
      POINTER (IZALFAT    , ZALFAT     (1) )
      POINTER (IZDLAT     , ZDLAT      (1) )
      POINTER (IZFCOR     , ZFCOR      (1) )
      POINTER (IZFTEMP    , ZFTEMP     (1) )
      POINTER (IZFVAP     , ZFVAP      (1) )
      POINTER (IZQDIAG    , ZQDIAG     (1) )
      POINTER (IZTDIAG    , ZTDIAG     (1) )
      POINTER (IZTSURF    , ZTSURF     (1) )
      POINTER (IZTSRAD    , ZTSRAD     (1) )
      POINTER (IZUDIAG    , ZUDIAG     (1) )
      POINTER (IZVDIAG    , ZVDIAG     (1) )
      POINTER (IZUE2      , ZUE2       (1) )
      POINTER (IZZA       , ZZA        (1) )
CZWQ
      POINTER (IZWW3LAT   , ZWW3LAT    (1) )
      POINTER (IZWW3LON   , ZWW3LON    (1) )
*
*
***
*
#include "zuzt.cdk"
*
      REAL Z0MAX, Z0HCON
#include "vamin.cdk"
      SAVE VAMIN, Z0MAX, Z0HCON
CZWQ
      logical seaspray,charnockl
*
***
*
***  WARNING ------ the value for Z0MAX needs to be increased when coupled with WAM
*
**    DATA Z0MAX / 5.E-2 /
**    DATA Z0MAX / 5.E-3 /
      DATA Z0HCON/ 4.E-5/
*
*** ------------------------------------------------------------------
*
#include "locbus.cdk"
      INTEGER INDX_SFC, SURFLEN
      PARAMETER (INDX_SFC = INDX_WATER)
      INTEGER QUELNIVO(MAXVARSURF)
*
#include "dintern.cdk"
*
*
#include "sfcbus.cdk"
*
#include "xptsurf.cdk"
*
#include "fintern.cdk"
*
      SURFLEN = M
*
*
*     EQUIVALENCES
*
      INIT_LOCBUS()
*
*     Syntax of macro locbus (must be typed in CAPITAL letters):
*     locbus (pointer, array_name_in_the_bus, level)
*     If level=0, array chosen automatically as follows:
*        1) level =  1 if array has  1 level only (e.g. TSURF )
*        2) level = nk if array has nk levels     (e.g. TMOINS)
*        3) level = indx_sfc if array has a level for each surface type (e.g. FC)
*        4) level has to be specified by user if array has more than one level
*           that all "belong" to the same surface type (e.g. TSOIL)
*
      LOCBUS (IALVIS_WAT, ALVIS  ,  0 )
      LOCBUS (ICMU      , BM     ,  0 )
      LOCBUS (ICTU      , BT     ,  0 )
      LOCBUS (IFC_WAT   , FC     ,  0 )
      LOCBUS (IFV_WAT   , FV     ,  0 )
      LOCBUS (IHST_WAT  , HST    ,  0 )
      LOCBUS (IHU       , HUMOINS,  0 )
      LOCBUS (IILMO_WAT , ILMO   ,  0 )
      LOCBUS (IPS       , PMOINS ,  0 )
      LOCBUS (IQS       , QSURF  ,  0 )
      LOCBUS (ITH       , THETAA ,  0 )
      LOCBUS (ITS       , TWATER ,  0 )
      LOCBUS (ITT       , TMOINS ,  0 )
      LOCBUS (IUU       , UMOINS ,  0 )
      LOCBUS (IVV       , VMOINS ,  0 )
      LOCBUS (IZ0H      , Z0T    ,  0 )
      LOCBUS (IZ0M      , Z0     ,  0 )
      LOCBUS (IZALFAQ   , ALFAQ  ,  0 )
      LOCBUS (IZALFAT   , ALFAT  ,  0 )
      LOCBUS (IZDLAT    , DLAT   ,  0 )
      LOCBUS (IZFCOR    , FCOR   ,  0 )
      LOCBUS (IZFTEMP   , FTEMP  ,  0 )
      LOCBUS (IZFVAP    , FVAP   ,  0 )
      LOCBUS (IZTSURF   , TSURF  ,  0 )
      LOCBUS (IZTSRAD   , TSRAD  ,  0 )
      LOCBUS (IZUDIAG   , UDIAG  ,  0 )
      LOCBUS (IZVDIAG   , VDIAG  ,  0 )
      LOCBUS (IZTDIAG   , TDIAG  ,  0 )
      LOCBUS (IZQDIAG   , QDIAG  ,  0 )
      LOCBUS (IZUE2     , UE2    ,  0 )
      LOCBUS (IZZA      , ZA     ,  0 )
CZWQ
      LOCBUS (IZWW3LAT  , WW3LAT ,  0 )
      LOCBUS (IZWW3LON  , WW3LON ,  0 )
*
      do i=1,n
         zun(i) = zu
         ztn(i) = zt
      end do
*
*------------------------------------------------------------------------
*
*
*       1.     Saturated specific humidity at the water surface
*       -------------------------------------------------------
*
*                           Uses FOQSA instead of FOQST to take into account saturation
*                           with respect to sea water (liquid between 0 and -1.8 C)
      DO I=1,N
        QS(I) = FOQSA(TS(I),PS(I))
      END DO
*
*
*
*       2.     Calculate roughness lengths based on generalized Charnock's relation 
*       ---------------------------------------------------------------------------
*
*
      beta = 0.018
*
       DO I=1,N
        VMOD  (I) = SQRT(MAX(VAMIN,UU(I)**2+VV(I)**2))
       END DO
*
       if (charnockl)  then
*        Z0MAX = 5.E-3 
        Z0MAX = 1.E-2 
        IF ( kount.gt.0 ) THEN
          DO I=1,N
           Z0M(I) = MAX( MIN( BETA*ZUE2(I)/GRAV,Z0MAX ) , Z0MIN )
          END DO
        ENDIF
       else
        Z0MAX = 1.E-2 
        IF (kount.gt.0 .and. kount .le. 3) THEN
           DO I=1,N
            Z0M(I) = MAX( MIN( BETA*ZUE2(I)/GRAV,Z0MAX ),Z0MIN )
           END DO
        ELSEIF (kount .gt. 3) THEN
           open ( 78,file='../cp.dat',form='unformatted')
           read(78)CPC
           read(78)WAGE
           close(78)
	      DO I=1,N
                 IC = int((ZWW3LON(I) - LONGS)/DX) +1
                 JC = int((ZWW3LAT(I) - LATS) /DY) +1
	         CP(I)   = CPC(IC,JC)
	         CHW3(I) = WAGE(IC,JC)
	       ENDDO
          CALL SMITH1(Z0M,CP,CHW3,ZUE2,VMOD,BETA,
     $                Z0MAX,Z0MIN,10,GRAV,N,.false.)     !true for wage from ww3
        ENDIF
      endif
*
*                           Arbitrary values for the height
*                           of the boundary layer (for FLXSURF1)
*                           No impact on the results.
      DO I=1,N
        IF (Z0TCST) THEN
           Z0H(I) = Z0HCON
        ELSE
           Z0H(I) = Z0M(I)
        ENDIF 
CZWQ        VMOD  (I) = SQRT(MAX(VAMIN,UU(I)**2+VV(I)**2))
      END DO
*
*
*
*       3.     Calculate the surface transfer coefficient and fluxes 
*       ------------------------------------------------------------
*
      IF ( .not. seaspray ) THEN 
      CALL FLXSURF2( CMU, CTU, SCR1, ZFTEMP, ZFVAP,
     $               ILMO_WAT, ZUE2, ZFCOR, TH, HU,
     $               ZZA, VMOD, TS, QS, HST_WAT,
     $               Z0M, Z0H,SCR2, SCR3, SCR4, 
     $               SCR5, SCR6, SCR7, N     )
      ELSE
      CALL FLXSURF3( CMU, CTU, SCR1, ZFTEMP, ZFVAP,
     $               ILMO_WAT, ZUE2, ZFCOR, TH, HU,
     $               ZZA, VMOD, TS, QS, HST_WAT,
     $               Z0M, Z0H,SCR2, SCR3, SCR4, 
     $               SCR5, SCR6, SCR7, N ,PS )
      ENDIF	 
*
*
      CALL DIASURF1(ZUDIAG, ZVDIAG, ZTDIAG, ZQDIAG,
     $              N, UU, VV, TS, QS,
     $              Z0M, Z0H, ILMO_WAT, ZZA,
     $              HST_WAT, ZUE2, ZFTEMP, ZFVAP,
     $              ZUN, ZTN, ZDLAT)
*
*
*
*       4.     Finalize the fluxes
*       --------------------------
*
*VDIR NODEP
      DO I=1,N
*
        ZTSURF   (I) = TS (I)
        ZTSRAD   (I) = TS (I)
*
****    calculs valides dans le cas explicite seulement
      IF ( .not. seaspray ) THEN 
        ZALFAT   (I) = - CTU(I) * ( TS(I)-TH(I) )
        ZALFAQ   (I) = - CTU(I) * ( QS(I)-HU(I) )
      ELSE 
        ZALFAT   (I) = -ZFTEMP(I)
        ZALFAQ   (I) = -ZFVAP(I)
      ENDIF
        CTU      (I) = 0.
        RHO = PS(I)/(RGASD * ZTDIAG(I)*(1.+DELTA*ZQDIAG(I)))
        FC   _WAT(I) = -CPD *RHO*ZALFAT(I)
        FV   _WAT(I) = -CHLC*RHO*ZALFAQ(I)
****
*
      END DO
*
*     FILL THE ARRAYS TO BE AGGREGATED LATER IN S/R AGREGE
      CALL FILLAGG ( BUS, BUSSIZ, PTSURF, PTSURFSIZ, INDX_WATER,
     +               SURFLEN )
*
      RETURN
      END
