ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90

MODULE sbcblk_mfs
   !!======================================================================
   !!                       ***  MODULE  sbcblk_mfs  ***
   !! Ocean forcing:  momentum, heat and freshwater flux formulation
   !!=====================================================================
   !! History :  3.3  !   2010-05 (P. Oddo) Original Code
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   sbc_blk_mfs  : bulk formulation as ocean surface boundary condition
   !!                   (forced mode, mfs bulk formulae)
   !!   blk_oce_mfs  : ocean: computes momentum, heat and freshwater fluxes
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE phycst          ! physical constants
   USE fldread         ! read input fields
   USE sbc_oce         ! Surface boundary condition: ocean fields
   USE iom             ! I/O manager library
   USE in_out_manager  ! I/O manager
   USE lib_mpp         ! distribued memory computing library
   USE wrk_nemo        ! work arrays
   USE timing          ! Timing
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
   USE prtctl          ! Print control
   USE sbcwave,ONLY : cdn_wave !wave module

   IMPLICIT NONE
   PRIVATE

   PUBLIC   sbc_blk_mfs       ! routine called in sbcmod module
      
   INTEGER , PARAMETER ::   jpfld   = 7         ! maximum number of files to read 
   INTEGER , PARAMETER ::   jp_wndi = 1         ! index of 10m wind velocity (i-component) (m/s) at T-point
   INTEGER , PARAMETER ::   jp_wndj = 2         ! index of 10m wind velocity (j-component) (m/s) at T-point
   INTEGER , PARAMETER ::   jp_clc  = 3         ! index of total cloud cover               ( % )
   INTEGER , PARAMETER ::   jp_msl  = 4         ! index of mean sea level pressure         (Pa)
   INTEGER , PARAMETER ::   jp_tair = 5         ! index of 10m air temperature             (Kelvin)
   INTEGER , PARAMETER ::   jp_rhm  = 6         ! index of dew point temperature           (Kelvin)
   INTEGER , PARAMETER ::   jp_prec = 7         ! index of total precipitation (rain+snow) (Kg/m2/s)
   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf ! structure of input fields (file informations, fields read)
         
   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
   !! $Id: sbcblk_mfs.F90 4578 2017-09-25 09:34:12Z ufla $
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS


   SUBROUTINE sbc_blk_mfs( kt )
      !!---------------------------------------------------------------------
      !!                    ***  ROUTINE sbc_blk_mfs  ***
      !!                   
      !! ** Purpose :   provide at each time step the surface ocean fluxes
      !!      (momentum, heat, freshwater, runoff is added later in the code) 
      !!
      !! ** Method  : (1) READ Atmospheric data from NetCDF files:
      !!      the 10m wind velocity (i-component) (m/s) at T-point
      !!      the 10m wind velocity (j-component) (m/s) at T-point
      !!      the 2m Dew point Temperature        (k)
      !!      the Cloud COver                     (%)
      !!      the 2m air temperature              (Kelvin)
      !!      the Mean Sea Level Preesure         (hPa)
      !!      the Climatological Precipitation    (kg/m2/s)
      !!              (2) CALL blk_oce_mfs
      !!
      !!      Computes:
      !!      Solar Radiation using Reed formula (1975, 1977)
      !!      Net Long wave radiation using Bignami et al. (1995)
      !!      Latent and Sensible heat using Kondo (1975)
      !!      Drag coeff using Hllerman and Rosenstein (1983)
      !!      C A U T I O N : never mask the surface stress fields
      !!                      the stress is assumed to be in the mesh referential
      !!                      i.e. the (i,j) referential
      !!
      !! ** Action  :   defined at each time-step at the air-sea interface
      !!              - utau, vtau  i- and j-component of the wind stress
      !!              - taum        wind stress module at T-point
      !!              - wndm        10m wind module at T-point over free ocean or leads in presence of sea-ice
      !!              - qns, qsr    non-slor and solar heat flux
      !!              - emp         evaporation minus precipitation
      !!----------------------------------------------------------------------
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sh_now   ! specific humidity at T-point 
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  catm     ! Cover 
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  alonl    ! Lon 
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  alatl    ! Lat 
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gsst     ! SST 
      !!---------------------------------------------------------------------
      !! Local fluxes variables
      !!---------------------------------------------------------------------
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  qbw     ! Net Long wave 
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ha      ! Sesnible 
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  elat    ! Latent 
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  evap    ! evaporation rate 

      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
      !!
      INTEGER  :: ierror                          ! return error code
      INTEGER  :: ifpr     ! dummy loop indice
      INTEGER  :: jj,ji    ! dummy loop arguments
      INTEGER  ::   ios    ! Local integer output status for namelist read
      REAL(wp) :: act_hour
      !!--------------------------------------------------------------------
      !! Variables for specific humidity computation
      !!--------------------------------------------------------------------
      REAL(wp) :: onsea,par1,par2
      DATA onsea,par1,par2 / 0.98, 640380., -5107.4 /
      !!                      par1 [Kg/m3], par2 [K]

      CHARACTER(len=100) ::  cn_dir                           ! Root directory for location of Atmospheric forcing files
      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                ! array of namelist informations on the fields to read
      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_clc, sn_msl       ! informations about the fields to be read
      TYPE(FLD_N) ::   sn_tair , sn_rhm, sn_prec              !   "                                 "
      !!---------------------------------------------------------------------
      NAMELIST/namsbc_mfs/ cn_dir ,                                          &
         &                  sn_wndi , sn_wndj, sn_clc   , sn_msl ,           &
         &                  sn_tair , sn_rhm , sn_prec 
      !!---------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('sbc_blk_mfs')
      !
      !                                         ! ====================== !
      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
         !                                      ! ====================== !
         ALLOCATE( sh_now(jpi,jpj), catm(jpi,jpj), alonl(jpi,jpj), alatl(jpi,jpj),     &
         &        gsst(jpi,jpj),  qbw(jpi,jpj),    ha(jpi,jpj),  elat(jpi,jpj),     &
         &        evap(jpi,jpj), STAT=ierror )

         IF( ierror /= 0 )   CALL ctl_warn('sbc_blk_mfs: failed to allocate arrays')

         REWIND( numnam_ref )              ! Namelist namsbc_msf in reference namelist : MFS files
         READ  ( numnam_ref, namsbc_mfs, IOSTAT = ios, ERR = 901)
901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_mfs in reference namelist', lwp )

         REWIND( numnam_cfg )              ! Namelist namsbc_msf in configuration namelist : MFS files
         READ  ( numnam_cfg, namsbc_mfs, IOSTAT = ios, ERR = 902 )
902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_mfs in configuration namelist', lwp )
         IF(lwm) WRITE ( numond, namsbc_mfs )
         !
         ! store namelist information in an array
         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
         slf_i(jp_clc ) = sn_clc    ;   slf_i(jp_msl ) = sn_msl
         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_rhm)  = sn_rhm
         slf_i(jp_prec) = sn_prec   ;  
         !
         ALLOCATE( sf(jpfld), STAT=ierror )         ! set sf structure
         IF( ierror > 0 ) THEN
            CALL ctl_stop( 'sbc_blk_mfs: unable to allocate sf structure' )   ;   RETURN
         ENDIF
         DO ifpr= 1, jpfld
            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
         END DO
         ! fill sf with slf_i and control print
         CALL fld_fill( sf, slf_i, cn_dir,'sbc_blk_mfs','bulk formulation for ocean SBC', 'namsbc_mfs' )
            !
      ENDIF

         CALL fld_read( kt, nn_fsbc, sf )                   ! input fields provided at the current time-step

         catm(:,:)   = 0.0    ! initializze cloud cover variable
         sh_now(:,:) = 0.0    ! initializze specifif humidity variable

         DO jj = 1, jpj
            DO ji = 1, jpi  

         ! Calculate Specific Humidity 
         !-------------------------------------------------
            sh_now(ji,jj) = (1/1.22) * onsea * par1 * EXP(par2/sf(jp_rhm)%fnow(ji,jj,1))

         ! Normalize Clouds
         !-------------------------------------------------
            catm(ji,jj)   = max(0.0,min(1.0,sf(jp_clc)%fnow(ji,jj,1)*0.01))

            END DO
         END DO

         ! wind module at 10m 
         !--------------------------------------------------
         wndm(:,:) = SQRT(  sf(jp_wndi)%fnow(:,:,1) * sf(jp_wndi)%fnow(:,:,1)   &
              &             + sf(jp_wndj)%fnow(:,:,1) * sf(jp_wndj)%fnow(:,:,1)  )

         ! Some conv for fluxes computation
         !-------------------------------------------------
         alonl(:,:) = glamt(:,:) * rad
         alatl(:,:) = gphit(:,:) * rad
         gsst(:,:)  = tsn(:,:,1,jp_tem)  * tmask(:,:,1)

         IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN

         ! Force to zero the output of fluxes 
         !-------------------------------------------------
          qsr(:,:)  = 0.0 ; qbw(:,:)  = 0.0 ; 
          ha(:,:)   = 0.0 ; elat(:,:) = 0.0 ; 
          evap(:,:) = 0.0 ; utau(:,:) = 0.0 ; 
          vtau(:,:) = 0.0

          CALL lbc_lnk( sf(jp_wndi)%fnow(:,:,1), 'T', -1. )
          CALL lbc_lnk( sf(jp_wndj)%fnow(:,:,1), 'T', -1. )

          act_hour = (( nsec_year / rday ) - INT (nsec_year / rday)) * rjjhh

          CALL fluxes_mfs(alatl,alonl,act_hour,                                &     ! input static
                            gsst(:,:),sf(jp_tair)%fnow(:,:,1),sh_now(:,:),     &     ! input dynamic
                            sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),  &     ! input dynamic
                            sf(jp_msl)%fnow(:,:,1) , catm(:,:) ,               &     ! input dynamic
                            qsr,qbw,ha,elat,evap,utau,vtau)                          ! output

         ! Shortwave radiation
         !--------------------------------------------------
          qsr(:,:) = qsr(:,:) * tmask(:,:,1)

         ! total non solar heat flux over water
         !--------------------------------------------------
          qns(:,:) = -1. * ( qbw(:,:) + ha(:,:) + elat(:,:) )
          qns(:,:) = qns(:,:)*tmask(:,:,1)

         ! mask the wind module at 10m 
         !--------------------------------------------------
          wndm(:,:) = wndm(:,:) * tmask(:,:,1)

         !   wind stress module (taum) into T-grid
         !--------------------------------------------------
          taum(:,:) = SQRT( utau(:,:) * utau(:,:) + vtau(:,:) * vtau(:,:) ) * tmask(:,:,1)

          CALL lbc_lnk( taum, 'T', 1. )

         ! Interpolate utau, vtau into the grid_V and grid_V
         !-------------------------------------------------
      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
         DO jj = 1, jpjm1
            DO ji = 1, fs_jpim1
               utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( utau(ji,jj) * tmask(ji,jj,1) &
               &                                + utau(ji+1,jj) * tmask(ji+1,jj,1) )        &
               &                 * MAX(tmask(ji,jj,1),tmask(ji+1,jj  ,1))
               vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vtau(ji,jj) * tmask(ji,jj,1) &
               &                                + vtau(ji,jj+1) * tmask(ji,jj+1,1) )        &
               &                 * MAX(tmask(ji,jj,1),tmask(ji  ,jj+1,1))
            END DO
         END DO

         CALL lbc_lnk( utau(:,:), 'U', -1. )
         CALL lbc_lnk( vtau(:,:), 'V', -1. )

         ! for basin budget and cooerence
         !--------------------------------------------------
!CDIR COLLAPSE
           emp (:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
!CDIR COLLAPSE

         CALL iom_put( "qlw_oce",   qbw  )                 ! output downward longwave heat over the ocean
         CALL iom_put( "qsb_oce", - ha   )                 ! output downward sensible heat over the ocean
         CALL iom_put( "qla_oce", - elat )                 ! output downward latent   heat over the ocean
         CALL iom_put( "qns_oce",   qns  )                 ! output downward non solar heat over the ocean

      ENDIF
      !
      IF( nn_timing == 1 )  CALL timing_stop('sbc_blk_mfs')
      !
   END SUBROUTINE sbc_blk_mfs
  
 
   SUBROUTINE fluxes_mfs(alat,alon,hour,                               &
        sst,tnow,shnow,unow,vnow,mslnow,cldnow,qsw,qbw,ha,elat,        &
        evap,taux,tauy)
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE fluxes_mfs  ***
      !!
      !! --- it provides SURFACE HEAT and MOMENTUM FLUXES in MKS :
      !!
      !!  1)   Water flux (WFLUX)                 [ watt/m*m ]
      !!  2)   Short wave flux (QSW)              [ watt/m*m ] Reed 1977
      !!  3)   Long wave flux backward (QBW)      [ watt/m*m ]
      !!  4)   Latent heat of evaporation (ELAT)  [ watt/m*m ]
      !!  5)   Sensible heat flux   (HA)          [ watt/m*m ]
      !!  6)   Wind stress x-component   (TAUX)   [ newton/m*m ]
      !!  7)   Wind stress y-component   (TAUY)   [ newton/m*m ]
      !!
      !!----------------------------------------------------------------------
      USE sbcblk_core, ONLY: turb_core_2z ! For wave coupling and Tair/rh from 2 to 10m

      REAL(wp), INTENT(in   ) :: hour
      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: sst, unow, alat , alon
      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: vnow, cldnow, mslnow
      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: qsw, qbw, ha, elat
      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: evap,taux,tauy
      REAL(wp), INTENT(inout), DIMENSION (:,:) :: tnow , shnow

      INTEGER :: ji,jj 
      REAL(wp)  :: wair, vtnow, ea, deltemp, s, stp , fh , fe
      REAL(wp)  :: esre, cseep

      REAL(wp), DIMENSION (:,:), POINTER ::   rspeed, sh10now, t10now, cdx, ce, shms
      REAL(wp), DIMENSION (:,:), POINTER ::   rhom, sstk, ch, rel_windu, rel_windv
      !!----------------------------------------------------------------------
      !!     coefficients ( in MKS )  :
      !!----------------------------------------------------------------------

      REAL(wp), PARAMETER ::  ps = 1013.    ! --- surface air pressure
      REAL(wp), PARAMETER ::  expsi=0.622   ! --- expsi
      REAL(wp), PARAMETER ::  rd=287.       ! --- dry air gas constant
      REAL(wp), PARAMETER ::  cp=1005.      ! --- specific heat capacity
      REAL(wp), PARAMETER ::  onsea=0.98    ! --- specific humidity factors
      REAL(wp), PARAMETER ::  par1=640380.  ! [Kg/m3]
      REAL(wp), PARAMETER ::  par2=-5107.4  ! [K]

      !---------------------------------------------------------------------
      !--- define Kondo parameters
      !---------------------------------------------------------------------

      REAL(wp), DIMENSION(5) :: a_h = (/0.0,0.927,1.15,1.17,1.652/)
      REAL(wp), DIMENSION(5) :: a_e = (/0.0,0.969,1.18,1.196,1.68/)
      REAL(wp), DIMENSION(5) :: b_h = (/1.185,0.0546,0.01,0.0075,-0.017/)
      REAL(wp), DIMENSION(5) :: b_e = (/1.23,0.0521,0.01,0.008,-0.016/)
      REAL(wp), DIMENSION(5) :: c_h = (/0.0,0.0,0.0,-0.00045,0.0/)
      REAL(wp), DIMENSION(5) :: c_e = (/0.0,0.0,0.0,-0.0004,0.0/)
      REAL(wp), DIMENSION(5) :: p_h = (/-0.157,1.0,1.0,1.0,1.0/)
      REAL(wp), DIMENSION(5) :: p_e = (/-0.16,1.0,1.0,1.0,1.0/)
      INTEGER :: kku                        !index varing with wind speed
      !
      IF( nn_timing == 1 )  CALL timing_start('fluxes_mfs')
      !
      CALL wrk_alloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
      CALL wrk_alloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )

      !!----------------------------------------------------------------------
      !! ------------------ (i)      short wave
      !!----------------------------------------------------------------------

       CALL qshort(hour,alat,alon,cldnow,qsw)

          rel_windu(:,:) = 0.0_wp
          rel_windv(:,:) = 0.0_wp

       DO jj = 2, jpj
          DO ji = fs_2, jpi
           rel_windu(ji,jj) = unow(ji,jj) - 0.5_wp * ( ssu_m(ji-1,jj) + ssu_m(ji,jj) )
           rel_windv(ji,jj) = vnow(ji,jj) - 0.5_wp * ( ssv_m(ji,jj-1) + ssv_m(ji,jj) )
          END DO
       END DO
     
       CALL lbc_lnk( rel_windu(:,:), 'U', -1. )
       CALL lbc_lnk( rel_windv(:,:), 'V', -1. )

       rspeed(:,:)= SQRT(rel_windu(:,:)*rel_windu(:,:)   &
         &                   + rel_windv(:,:)*rel_windv(:,:)) 

       sstk(:,:) = sst(:,:) + rtt                          !- SST data in Kelvin degrees
       shms(:,:) = (1/1.22)*onsea*par1*EXP(par2/sstk(:,:)) !- Saturation Specific Humidity

      ! --- Transport temperature and humidity from 2m to 10m
      !----------------------------------------------------------------------

      t10now(:,:) = 0.0           ;   sh10now(:,:)= 0.0
      ! Note that air temp is converted in potential temp
      CALL turb_core_2z(2.,10.,sstk,tnow+2*0.0098,shms,shnow,rspeed,        &
         &              Cdx,Ch,Ce,t10now,sh10now )
      tnow(:,:)  = t10now(:,:)    ;    shnow(:,:) = sh10now(:,:)

      !!----------------------------------------------------------------------
      !! ------------------ (ii)    net long wave
      !!----------------------------------------------------------------------

      DO jj = 1, jpj
         DO ji = 1, jpi
            wair = shnow(ji,jj) / (1 - shnow(ji,jj))    ! mixing ratio of the air (Wallace and Hobbs)
            vtnow = (tnow(ji,jj)*(expsi+wair))/(expsi*(1.+wair))   ! virtual temperature of air
            rhom(ji,jj) = 100.*(ps/rd)/vtnow                       ! density of the moist air
            ea   = (wair  / (wair  + 0.622 )) * mslnow(ji,jj)

            qbw(ji,jj) = emic*stefan*( sstk(ji,jj)**4. )                    &
                 - ( stefan*( tnow(ji,jj)**4. ) * ( 0.653 + 0.00535*ea ) )  &
                   * ( 1. + 0.1762*( cldnow(ji,jj)**2. ) )

         END DO
      END DO

      DO jj = 1, jpj
         DO ji = 1, jpi
      !!----------------------------------------------------------------------
      !! ------------------ (iii)   sensible heat
      !!----------------------------------------------------------------------

      !! --- calculates the term :      ( Ts - Ta )
      !!----------------------------------------------------------------------
            deltemp = sstk(ji,jj) - tnow (ji,jj)

      !!----------------------------------------------------------------------
      !! --- variable turbulent exchange coefficients ( from Kondo 1975 )
      !! --- calculate the Neutral Transfer Coefficent using an empiric formula
      !! --- by Kondo et al. Then it applies the diabatic approximation.
      !!----------------------------------------------------------------------

            s = deltemp/(wndm(ji,jj)**2.)   !! --- calculate S 
            stp = s*abs(s)/(abs(s)+0.01)    !! --- calculate the Stability Parameter 

      !!----------------------------------------------------------------------
      !! --- for stable condition (sst-t_air < 0):
      !!----------------------------------------------------------------------

            IF (s.lt.0. .and. ((stp.gt.-3.3).and.(stp.lt.0.))) THEN
                fh = 0.1_wp+0.03_wp*stp+0.9_wp*exp(4.8_wp*stp)
                fe = fh
            ELSE IF (s.lt.0. .and. stp.le.-3.3) THEN
                fh = 0._wp
                fe = fh
            ELSE                                       ! --- for unstable condition 
                fh = 1.0_wp+0.63_wp*sqrt(stp)
                fe = fh
            ENDIF

      !!----------------------------------------------------------------------
      !! --- calculate the coefficient CH,CE,CD
      !!----------------------------------------------------------------------

            IF (wndm(ji,jj) >= 0. .AND. wndm(ji,jj) <= 2.2)       THEN
                kku=1
            ELSE IF (wndm(ji,jj) > 2.2 .AND. wndm(ji,jj) <= 5.0)  THEN
                kku=2
            ELSE IF (wndm(ji,jj) > 5.0 .AND. wndm(ji,jj) <= 8.0)  THEN
                kku=3
            ELSE IF (wndm(ji,jj) > 8.0 .AND. wndm(ji,jj) <= 25.0) THEN
                kku=4
            ELSE IF (wndm(ji,jj) > 25.0 )                         THEN
                kku=5
            ENDIF

            ch(ji,jj) = ( a_h(kku) + b_h(kku) * wndm(ji,jj) ** p_h(kku)      &
                        + c_h(kku) * (wndm(ji,jj)-8 ) **2) * fh

            ce(ji,jj) = ( a_e(kku) + b_e(kku) * wndm(ji,jj) ** p_e(kku)      &
                        + c_e(kku) * (wndm(ji,jj)-8 ) **2) * fe

            ch(ji,jj) = ch(ji,jj) / 1000.0
            ce(ji,jj) = ce(ji,jj) / 1000.0

            IF (wndm(ji,jj)<0.3) THEN
               ch(ji,jj) = 1.3e-03 * fh
               ce(ji,jj) = 1.5e-03 * fe
            ELSE IF(wndm(ji,jj)>50.0) THEN
               ch(ji,jj) = 1.25e-03 * fh
               ce(ji,jj) = 1.30e-03 * fe
            ENDIF

      !!----------------------------------------------------------------------
      !! --- calculates the SENSIBLE HEAT FLUX in MKS ( watt/m*m )
      !!----------------------------------------------------------------------

            HA(ji,jj) = rhom(ji,jj)*cp*ch(ji,jj)*wndm(ji,jj)*deltemp

      !!----------------------------------------------------------------------
      !! ------------------ (iv)  latent heat
      !! --- calculates the LATENT HEAT FLUX  ( watt/m*m )
      !! --- ELAT = L*rho*Ce*|V|*[qs(Ts)-qa(t2d)]
      !!----------------------------------------------------------------------

            esre  = shms(ji,jj) - shnow(ji,jj)   ! --- calculates the term : qs(Ta)-qa(t2d)

            cseep = ce(ji,jj) * wndm(ji,jj) * esre     ! --- calculates the term : Ce*|V|*[qs(Ts)-qa(t2d)]

            evap(ji,jj) = (cseep * rhom(ji,jj))  ! in [kg/m2/sec] !! --- calculates the EVAPORATION RATE [m/yr]

            elat(ji,jj) = rhom(ji,jj) * cseep * heatlat(sst(ji,jj))

      !!----------------------------------------------------------------------
      !! --- calculates the Drag Coefficient
      !!----------------------------------------------------------------------

      !!----------------------------------------------------------------------
      !! --- deltemp should be (Ts - Ta) in the formula estimating
      !! --- drag coefficient
      !!----------------------------------------------------------------------

              IF( .NOT. ln_cdgw ) THEN
                 cdx(ji,jj) = cd_HR(wndm(ji,jj),deltemp)
              ENDIF

          END DO
      END DO

      !!----------------------------------------------------------------------
      !! --- calculates the wind stresses in MKS ( newton/m*m )
      !! ---            taux= rho*Cd*|V|u      tauy= rho*Cd*|V|v
      !!----------------------------------------------------------------------

       taux(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windu(:,:)
       tauy(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windv(:,:)

      CALL wrk_dealloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
      CALL wrk_dealloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
      !
      IF( nn_timing == 1 )  CALL timing_stop('fluxes_mfs')
      !
   END SUBROUTINE fluxes_mfs


      REAL(wp) FUNCTION cd_HR(speed,delt)
      !!----------------------------------------------------------------------
      !! --- calculates the Drag Coefficient as a function of the abs. value
      !! --- of the wind velocity ( Hellermann and Rosenstein )
      !!----------------------------------------------------------------------

       REAL(wp), INTENT(in) :: speed,delt
       REAL(wp), PARAMETER  :: a1=0.934e-3 , a2=0.788e-4, a3=0.868e-4     
       REAL(wp), PARAMETER  :: a4=-0.616e-6, a5=-.120e-5, a6=-.214e-5

        cd_HR = a1 + a2*speed + a3*delt + a4*speed*speed        &
           + a5*delt*delt  + a6*speed*delt

      END FUNCTION cd_HR

      REAL(wp) function HEATLAT(t)
      !!----------------------------------------------------------------------
      !! --- calculates the Latent Heat of Vaporization ( J/kg ) as function of
      !! --- the temperature ( Celsius degrees )
      !! --- ( from A. Gill  pag. 607 )
      !!
      !! --- Constant Latent Heat of Vaporization ( Rosati,Miyakoda 1988 )
      !!     L = 2.501e+6  (MKS)
      !!----------------------------------------------------------------------

      REAL(wp) , intent(in) :: t

      heatlat = 2.5008e+6 -2.3e+3*t

      END FUNCTION HEATLAT


   SUBROUTINE qshort(hour,alat,alon,cldnow,qsw)
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE qshort  ***
      !!
      !! ** Purpose :   Compute Solar Radiation
      !!
      !! ** Method  :   Compute Solar Radiation according Astronomical
      !!                formulae
      !!
      !! References :   Reed RK (1975) and Reed RK (1977)
      !!
      !! Note: alat,alon - (lat, lon)  in radians
      !!----------------------------------------------------------------------
        REAL(wp), INTENT (in) :: hour

        REAL(wp), INTENT(in ), DIMENSION(:,:) :: alat,alon
        REAL(wp), INTENT(in ), DIMENSION(:,:) :: cldnow
        REAL(wp), INTENT(out), DIMENSION(:,:) :: qsw
        REAL(wp), DIMENSION(12) :: alpham

        REAL(wp), PARAMETER ::   eclips=23.439* (3.141592653589793_wp / 180._wp)
        REAL(wp), PARAMETER ::   solar = 1350.
        REAL(wp), PARAMETER ::   tau = 0.7
        REAL(wp), PARAMETER ::   aozone = 0.09
        REAL(wp), PARAMETER ::   yrdays = 360.
        REAL(wp) :: days, th0,th02,th03, sundec, thsun, coszen, qatten
        REAL(wp) :: qzer, qdir,qdiff,qtot,tjul,sunbet
        REAL(wp) :: albedo
        INTEGER :: jj, ji

      !!----------------------------------------------------------------------
      !! --- albedo monthly values from Payne (1972) as means of the values
      !! --- at 40N and 30N for the Atlantic Ocean ( hence the same latitudinal
      !! --- band of the Mediterranean Sea ) :
      !!----------------------------------------------------------------------

        data alpham /0.095,0.08,0.065,0.065,0.06,0.06,0.06,0.06,        &
                    0.065,0.075,0.09,0.10/

      !!----------------------------------------------------------------------
      !!   days is the number of days elapsed until the day=nday_year
      !!----------------------------------------------------------------------
        days = nday_year -1.
        th0  = 2.*rpi*days/yrdays
        th02 = 2.*th0
        th03 = 3.*th0

      !! --- sun declination :
      !!----------------------------------------------------------------------
        sundec = 0.006918 - 0.399912*cos(th0) + 0.070257*sin(th0) -   &
                          0.006758*cos(th02) + 0.000907*sin(th02) -   &
                          0.002697*cos(th03) + 0.001480*sin(th03)

      DO jj = 1, jpj
         DO ji = 1, jpi

      !! --- sun hour angle :
      !!----------------------------------------------------------------------
          thsun = (hour -12.)*15.*rad + alon(ji,jj)

      !! --- cosine of the solar zenith angle :
      !!----------------------------------------------------------------------
          coszen =sin(alat(ji,jj))*sin(sundec)                 &
                    +cos(alat(ji,jj))*cos(sundec)*cos(thsun)

          IF(coszen .LE. 5.035D-04) THEN
            coszen = 0.0
            qatten = 0.0
          ELSE
            qatten = tau**(1./coszen)
          END IF

          qzer  = coszen * solar *                                 &
                  (1.0+1.67E-2*cos(rpi*2.*(days-3.0)/365.0))**2
          qdir  = qzer * qatten
          qdiff = ((1.-aozone)*qzer - qdir) * 0.5
          qtot  =  qdir + qdiff
          tjul = (days -81.)*rad

      !! --- sin of the solar noon altitude in radians :
      !!----------------------------------------------------------------------
          sunbet=sin(alat(ji,jj))*sin(eclips*sin(tjul)) +   &
                 cos(alat(ji,jj))*cos(eclips*sin(tjul))

      !! --- solar noon altitude in degrees :
      !!----------------------------------------------------------------------

         sunbet = asin(sunbet)/rad

      !!----------------------------------------------------------------------
      !! --- calculates the albedo according to Payne (1972)
      !!----------------------------------------------------------------------

         albedo = alpham(nmonth)

      !!----------------------------------------------------------------------
      !! --- ( radiation as from Reed(1977), Simpson and Paulson(1979) )
      !! --- calculates SHORT WAVE FLUX ( watt/m*m )
      !! --- ( Rosati,Miyakoda 1988 ; eq. 3.8 )
      !!----------------------------------------------------------------------

          IF(cldnow(ji,jj).LT.0.3) THEN
             qsw(ji,jj) = qtot * (1.-albedo)
          ELSE
             qsw(ji,jj) = qtot*(1.-0.62*cldnow(ji,jj)              &
                                + .0019*sunbet)*(1.-albedo)
          ENDIF

         END DO
      END DO

   END SUBROUTINE qshort


   !!======================================================================

END MODULE sbcblk_mfs