New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcblk_mfs.F90 in branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 7380

Last change on this file since 7380 was 7380, checked in by emanuelaclementi, 7 years ago

#1805 updated 2016/dev_INGV_UKMO_2016
  • Property svn:keywords set to Id
File size: 29.8 KB
Line 
1MODULE sbcblk_mfs
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk_mfs  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!=====================================================================
6   !! History :  3.3  !   2010-05 (P. Oddo) Original Code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   sbc_blk_mfs   : bulk formulation as ocean surface boundary condition
11   !!                   (forced mode, mfs bulk formulae)
12   !!   blk_oce_mfs   : ocean: computes momentum, heat and freshwater fluxes
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers
15   USE dom_oce        ! ocean space and time domain
16   USE phycst         ! physical constants
17   USE fldread        ! read input fields
18   USE sbc_oce        ! Surface boundary condition: ocean fields
19   !
20   USE iom            ! I/O manager library
21   USE in_out_manager ! I/O manager
22   USE lib_mpp        ! distribued memory computing library
23   USE wrk_nemo       ! work arrays
24   USE timing         ! Timing
25   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
26   USE prtctl         ! Print control
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   sbc_blk_mfs       ! routine called in sbcmod module
32     
33   INTEGER , PARAMETER ::   jpfld   = 7         ! maximum number of files to read
34   INTEGER , PARAMETER ::   jp_wndi = 1         ! index of 10m wind velocity (i-component) (m/s) at T-point
35   INTEGER , PARAMETER ::   jp_wndj = 2         ! index of 10m wind velocity (j-component) (m/s) at T-point
36   INTEGER , PARAMETER ::   jp_clc  = 3         ! index of total cloud cover               ( % )
37   INTEGER , PARAMETER ::   jp_msl  = 4         ! index of mean sea level pressure         (Pa)
38   INTEGER , PARAMETER ::   jp_tair = 5         ! index of 10m air temperature             (Kelvin)
39   INTEGER , PARAMETER ::   jp_rhm  = 6         ! index of dew point temperature           (Kelvin)
40   INTEGER , PARAMETER ::   jp_prec = 7         ! index of total precipitation (rain+snow) (Kg/m2/s)
41   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf ! structure of input fields (file informations, fields read)
42         
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
47   !! $Id$
48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE sbc_blk_mfs( kt )
53      !!---------------------------------------------------------------------
54      !!                    ***  ROUTINE sbc_blk_mfs  ***
55      !!                   
56      !! ** Purpose :   provide at each time step the surface ocean fluxes
57      !!      (momentum, heat, freshwater, runoff is added later in the code)
58      !!
59      !! ** Method  : (1) READ Atmospheric data from NetCDF files:
60      !!      the 10m wind velocity (i-component) (m/s) at T-point
61      !!      the 10m wind velocity (j-component) (m/s) at T-point
62      !!      the 2m Dew point Temperature        (k)
63      !!      the Cloud COver                     (%)
64      !!      the 2m air temperature              (Kelvin)
65      !!      the Mean Sea Level Preesure         (hPa)
66      !!      the Climatological Precipitation    (kg/m2/s)
67      !!              (2) CALL blk_oce_mfs
68      !!
69      !!      Computes:
70      !!      Solar Radiation using Reed formula (1975, 1977)
71      !!      Net Long wave radiation using Bignami et al. (1995)
72      !!      Latent and Sensible heat using Kondo (1975)
73      !!      Drag coeff using Hllerman and Rosenstein (1983)
74      !!      C A U T I O N : never mask the surface stress fields
75      !!                      the stress is assumed to be in the mesh referential
76      !!                      i.e. the (i,j) referential
77      !!
78      !! ** Action  :   defined at each time-step at the air-sea interface
79      !!              - utau, vtau  i- and j-component of the wind stress
80      !!              - taum        wind stress module at T-point
81      !!              - wndm        10m wind module at T-point over free ocean or leads in presence of sea-ice
82      !!              - qns, qsr    non-slor and solar heat flux
83      !!              - emp         evaporation minus precipitation
84      !!----------------------------------------------------------------------
85      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sh_now   ! specific humidity at T-point
86      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  catm     ! Cover
87      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  alonl    ! Lon
88      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  alatl    ! Lat
89      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gsst     ! SST
90      !!---------------------------------------------------------------------
91      !! Local fluxes variables
92      !!---------------------------------------------------------------------
93      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  qbw     ! Net Long wave
94      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ha      ! Sesnible
95      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  elat    ! Latent
96      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  evap    ! evaporation rate
97
98      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
99      !!
100      INTEGER  :: ierror                          ! return error code
101      INTEGER  :: ifpr     ! dummy loop indice
102      INTEGER  :: jj,ji    ! dummy loop arguments
103      INTEGER  ::   ios    ! Local integer output status for namelist read
104      REAL(wp) :: act_hour
105      !!--------------------------------------------------------------------
106      !! Variables for specific humidity computation
107      !!--------------------------------------------------------------------
108      REAL(wp) :: onsea,par1,par2
109      DATA onsea,par1,par2 / 0.98, 640380., -5107.4 /
110      !!                      par1 [Kg/m3], par2 [K]
111
112      CHARACTER(len=100) ::  cn_dir                           ! Root directory for location of Atmospheric forcing files
113      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                ! array of namelist informations on the fields to read
114      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_clc, sn_msl       ! informations about the fields to be read
115      TYPE(FLD_N) ::   sn_tair , sn_rhm, sn_prec              !   "                                 "
116      !!---------------------------------------------------------------------
117      NAMELIST/namsbc_mfs/ cn_dir ,                                          &
118         &                  sn_wndi , sn_wndj, sn_clc   , sn_msl ,           &
119         &                  sn_tair , sn_rhm , sn_prec 
120      !!---------------------------------------------------------------------
121      !
122      IF( nn_timing == 1 )  CALL timing_start('sbc_blk_mfs')
123      !
124      !                                         ! ====================== !
125      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
126         !                                      ! ====================== !
127         ALLOCATE( sh_now(jpi,jpj), catm(jpi,jpj), alonl(jpi,jpj), alatl(jpi,jpj),     &
128         &        gsst(jpi,jpj),  qbw(jpi,jpj),    ha(jpi,jpj),  elat(jpi,jpj),     &
129         &        evap(jpi,jpj), STAT=ierror )
130
131         IF( ierror /= 0 )   CALL ctl_warn('sbc_blk_mfs: failed to allocate arrays')
132
133         REWIND( numnam_ref )              ! Namelist namsbc_msf in reference namelist : MFS files
134         READ  ( numnam_ref, namsbc_mfs, IOSTAT = ios, ERR = 901)
135901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_mfs in reference namelist', lwp )
136
137         REWIND( numnam_cfg )              ! Namelist namsbc_msf in configuration namelist : MFS files
138         READ  ( numnam_cfg, namsbc_mfs, IOSTAT = ios, ERR = 902 )
139902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_mfs in configuration namelist', lwp )
140         IF(lwm) WRITE ( numond, namsbc_mfs )
141         !
142         ! store namelist information in an array
143         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
144         slf_i(jp_clc ) = sn_clc    ;   slf_i(jp_msl ) = sn_msl
145         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_rhm)  = sn_rhm
146         slf_i(jp_prec) = sn_prec   ; 
147         !
148         ALLOCATE( sf(jpfld), STAT=ierror )         ! set sf structure
149         IF( ierror > 0 ) THEN
150            CALL ctl_stop( 'sbc_blk_mfs: unable to allocate sf structure' )   ;   RETURN
151         ENDIF
152         DO ifpr= 1, jpfld
153            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
154            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
155         END DO
156         ! fill sf with slf_i and control print
157         CALL fld_fill( sf, slf_i, cn_dir,'sbc_blk_mfs','bulk formulation for ocean SBC', 'namsbc_mfs' )
158            !
159      ENDIF
160
161         CALL fld_read( kt, nn_fsbc, sf )                   ! input fields provided at the current time-step
162
163         catm(:,:)   = 0.0    ! initializze cloud cover variable
164         sh_now(:,:) = 0.0    ! initializze specifif humidity variable
165
166         DO jj = 2, jpjm1
167            DO ji = fs_2, fs_jpim1   ! vector opt.
168
169         ! Calculate Specific Humidity
170         !-------------------------------------------------
171            sh_now(ji,jj) = (1/1.22) * onsea * par1 * EXP(par2/sf(jp_rhm)%fnow(ji,jj,1))
172
173         ! Normalize Clouds
174         !-------------------------------------------------
175            catm(ji,jj)   = max(0.0,min(1.0,sf(jp_clc)%fnow(ji,jj,1)*0.01))
176
177            END DO
178         END DO
179
180         ! wind module at 10m
181         !--------------------------------------------------
182         wndm(:,:) = SQRT(  sf(jp_wndi)%fnow(:,:,1) * sf(jp_wndi)%fnow(:,:,1)   &
183              &             + sf(jp_wndj)%fnow(:,:,1) * sf(jp_wndj)%fnow(:,:,1)  )
184
185         ! Some conv for fluxes computation
186         !-------------------------------------------------
187         alonl(:,:) = glamt(:,:) * rad
188         alatl(:,:) = gphit(:,:) * rad
189         gsst(:,:)  = tsn(:,:,1,jp_tem)  * tmask(:,:,1)
190
191         IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
192
193         ! Force to zero the output of fluxes
194         !-------------------------------------------------
195          qsr(:,:)  = 0.0 ; qbw(:,:)  = 0.0 ; 
196          ha(:,:)   = 0.0 ; elat(:,:) = 0.0 ; 
197          evap(:,:) = 0.0 ; utau(:,:) = 0.0 ; 
198          vtau(:,:) = 0.0
199
200          CALL lbc_lnk( sf(jp_wndi)%fnow(:,:,1), 'T', -1. )
201          CALL lbc_lnk( sf(jp_wndj)%fnow(:,:,1), 'T', -1. )
202
203          act_hour = (( nsec_year / rday ) - INT (nsec_year / rday)) * rjjhh
204
205          CALL fluxes_mfs(alatl,alonl,act_hour,                                &     ! input static
206                            gsst(:,:),sf(jp_tair)%fnow(:,:,1),sh_now(:,:),     &     ! input dynamic
207                            sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),  &     ! input dynamic
208                            sf(jp_msl)%fnow(:,:,1) , catm(:,:) ,               &     ! input dynamic
209                            qsr,qbw,ha,elat,evap,utau,vtau)                          ! output
210
211         ! Shortwave radiation
212         !--------------------------------------------------
213          qsr(:,:) = qsr(:,:) * tmask(:,:,1)
214
215         ! total non solar heat flux over water
216         !--------------------------------------------------
217          qns(:,:) = -1. * ( qbw(:,:) + ha(:,:) + elat(:,:) )
218          qns(:,:) = qns(:,:)*tmask(:,:,1)
219
220         ! mask the wind module at 10m
221         !--------------------------------------------------
222          wndm(:,:) = wndm(:,:) * tmask(:,:,1)
223
224         !   wind stress module (taum) into T-grid
225         !--------------------------------------------------
226          taum(:,:) = SQRT( utau(:,:) * utau(:,:) + vtau(:,:) * vtau(:,:) ) * tmask(:,:,1)
227
228          CALL lbc_lnk( taum, 'T', 1. )
229
230         ! Interpolate utau, vtau into the grid_V and grid_V
231         !-------------------------------------------------
232      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
233      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
234         DO jj = 1, jpjm1
235            DO ji = 1, fs_jpim1
236               utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( utau(ji,jj) * tmask(ji,jj,1) &
237               &                                + utau(ji+1,jj) * tmask(ji+1,jj,1) )        &
238               &                 * MAX(tmask(ji,jj,1),tmask(ji+1,jj  ,1))
239               vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vtau(ji,jj) * tmask(ji,jj,1) &
240               &                                + vtau(ji,jj+1) * tmask(ji,jj+1,1) )        &
241               &                 * MAX(tmask(ji,jj,1),tmask(ji  ,jj+1,1))
242            END DO
243         END DO
244
245         CALL lbc_lnk( utau(:,:), 'U', -1. )
246         CALL lbc_lnk( vtau(:,:), 'V', -1. )
247
248         ! for basin budget and cooerence
249         !--------------------------------------------------
250         emp(:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
251
252         CALL iom_put( "qlw_oce",   qbw  )                 ! output downward longwave heat over the ocean
253         CALL iom_put( "qsb_oce", - ha   )                 ! output downward sensible heat over the ocean
254         CALL iom_put( "qla_oce", - elat )                 ! output downward latent   heat over the ocean
255         CALL iom_put( "qns_oce",   qns  )                 ! output downward non solar heat over the ocean
256
257      ENDIF
258      !
259      IF( nn_timing == 1 )  CALL timing_stop('sbc_blk_mfs')
260      !
261   END SUBROUTINE sbc_blk_mfs
262 
263 
264   SUBROUTINE fluxes_mfs(alat,alon,hour,                               &
265        sst,tnow,shnow,unow,vnow,mslnow,cldnow,qsw,qbw,ha,elat,        &
266        evap,taux,tauy)
267      !!----------------------------------------------------------------------
268      !!                    ***  ROUTINE fluxes_mfs  ***
269      !!
270      !! --- it provides SURFACE HEAT and MOMENTUM FLUXES in MKS :
271      !!
272      !!  1)   Water flux (WFLUX)                 [ watt/m*m ]
273      !!  2)   Short wave flux (QSW)              [ watt/m*m ] Reed 1977
274      !!  3)   Long wave flux backward (QBW)      [ watt/m*m ]
275      !!  4)   Latent heat of evaporation (ELAT)  [ watt/m*m ]
276      !!  5)   Sensible heat flux   (HA)          [ watt/m*m ]
277      !!  6)   Wind stress x-component   (TAUX)   [ newton/m*m ]
278      !!  7)   Wind stress y-component   (TAUY)   [ newton/m*m ]
279      !!
280      !!----------------------------------------------------------------------
281      USE sbcblk_core, ONLY: turb_core_2z ! For wave coupling and Tair/rh from 2 to 10m
282
283      REAL(wp), INTENT(in   ) :: hour
284      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: sst, unow, alat , alon
285      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: vnow, cldnow, mslnow
286      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: qsw, qbw, ha, elat
287      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: evap,taux,tauy
288      REAL(wp), INTENT(inout), DIMENSION (:,:) :: tnow , shnow
289
290      INTEGER :: ji,jj 
291      REAL(wp)  :: wair, vtnow, ea, deltemp, s, stp , fh , fe
292      REAL(wp)  :: esre, cseep
293
294      REAL(wp), DIMENSION (:,:), POINTER ::   rspeed, sh10now, t10now, cdx, ce, shms
295      REAL(wp), DIMENSION (:,:), POINTER ::   rhom, sstk, ch, rel_windu, rel_windv
296      !!----------------------------------------------------------------------
297      !!     coefficients ( in MKS )  :
298      !!----------------------------------------------------------------------
299
300      REAL(wp), PARAMETER ::  ps = 1013.    ! --- surface air pressure
301      REAL(wp), PARAMETER ::  expsi=0.622   ! --- expsi
302      REAL(wp), PARAMETER ::  rd=287.       ! --- dry air gas constant
303      REAL(wp), PARAMETER ::  cp=1005.      ! --- specific heat capacity
304      REAL(wp), PARAMETER ::  onsea=0.98    ! --- specific humidity factors
305      REAL(wp), PARAMETER ::  par1=640380.  ! [Kg/m3]
306      REAL(wp), PARAMETER ::  par2=-5107.4  ! [K]
307
308      !---------------------------------------------------------------------
309      !--- define Kondo parameters
310      !---------------------------------------------------------------------
311
312      REAL(wp), DIMENSION(5) :: a_h = (/0.0,0.927,1.15,1.17,1.652/)
313      REAL(wp), DIMENSION(5) :: a_e = (/0.0,0.969,1.18,1.196,1.68/)
314      REAL(wp), DIMENSION(5) :: b_h = (/1.185,0.0546,0.01,0.0075,-0.017/)
315      REAL(wp), DIMENSION(5) :: b_e = (/1.23,0.0521,0.01,0.008,-0.016/)
316      REAL(wp), DIMENSION(5) :: c_h = (/0.0,0.0,0.0,-0.00045,0.0/)
317      REAL(wp), DIMENSION(5) :: c_e = (/0.0,0.0,0.0,-0.0004,0.0/)
318      REAL(wp), DIMENSION(5) :: p_h = (/-0.157,1.0,1.0,1.0,1.0/)
319      REAL(wp), DIMENSION(5) :: p_e = (/-0.16,1.0,1.0,1.0,1.0/)
320      INTEGER :: kku                        !index varing with wind speed
321      !
322      IF( nn_timing == 1 )  CALL timing_start('fluxes_mfs')
323      !
324      CALL wrk_alloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
325      CALL wrk_alloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
326
327      !!----------------------------------------------------------------------
328      !! ------------------ (i)      short wave
329      !!----------------------------------------------------------------------
330
331       CALL qshort(hour,alat,alon,cldnow,qsw)
332
333          rel_windu(:,:) = 0.0_wp
334          rel_windv(:,:) = 0.0_wp
335
336       DO jj = 2, jpj
337          DO ji = fs_2, jpi
338           rel_windu(ji,jj) = unow(ji,jj) - 0.5_wp * ( ssu_m(ji-1,jj) + ssu_m(ji,jj) )
339           rel_windv(ji,jj) = vnow(ji,jj) - 0.5_wp * ( ssv_m(ji,jj-1) + ssv_m(ji,jj) )
340          END DO
341       END DO
342
343       rspeed(:,:)= SQRT(rel_windu(:,:)*rel_windu(:,:)   &
344         &                   + rel_windv(:,:)*rel_windv(:,:)) 
345
346       sstk(:,:) = sst(:,:) + rtt                          !- SST data in Kelvin degrees
347       shms(:,:) = (1/1.22)*onsea*par1*EXP(par2/sstk(:,:)) !- Saturation Specific Humidity
348
349      ! --- Transport temperature and humidity from 2m to 10m
350      !----------------------------------------------------------------------
351
352      t10now(:,:) = 0.0           ;   sh10now(:,:)= 0.0
353      ! Note that air temp is converted in potential temp
354      CALL turb_core_2z(2.,10.,sstk,tnow+2*0.0098,shms,shnow,rspeed,        &
355         &              Cdx,Ch,Ce,t10now,sh10now )
356      tnow(:,:)  = t10now(:,:)    ;    shnow(:,:) = sh10now(:,:)
357
358      !!----------------------------------------------------------------------
359      !! ------------------ (ii)    net long wave
360      !!----------------------------------------------------------------------
361
362      DO jj = 1, jpj
363         DO ji = 1, jpi
364            wair = shnow(ji,jj) / (1 - shnow(ji,jj))    ! mixing ratio of the air (Wallace and Hobbs)
365            vtnow = (tnow(ji,jj)*(expsi+wair))/(expsi*(1.+wair))   ! virtual temperature of air
366            rhom(ji,jj) = 100.*(ps/rd)/vtnow                       ! density of the moist air
367            ea   = (wair  / (wair  + 0.622 )) * mslnow(ji,jj)
368
369            qbw(ji,jj) = emic*stefan*( sstk(ji,jj)**4. )                    &
370                 - ( stefan*( tnow(ji,jj)**4. ) * ( 0.653 + 0.00535*ea ) )  &
371                   * ( 1. + 0.1762*( cldnow(ji,jj)**2. ) )
372
373         END DO
374      END DO
375
376      DO jj = 1, jpj
377         DO ji = 1, jpi
378      !!----------------------------------------------------------------------
379      !! ------------------ (iii)   sensible heat
380      !!----------------------------------------------------------------------
381
382      !! --- calculates the term :      ( Ts - Ta )
383      !!----------------------------------------------------------------------
384            deltemp = sstk(ji,jj) - tnow (ji,jj)
385
386      !!----------------------------------------------------------------------
387      !! --- variable turbulent exchange coefficients ( from Kondo 1975 )
388      !! --- calculate the Neutral Transfer Coefficent using an empiric formula
389      !! --- by Kondo et al. Then it applies the diabatic approximation.
390      !!----------------------------------------------------------------------
391
392            s = deltemp/(wndm(ji,jj)**2.)   !! --- calculate S
393            stp = s*abs(s)/(abs(s)+0.01)    !! --- calculate the Stability Parameter
394
395      !!----------------------------------------------------------------------
396      !! --- for stable condition (sst-t_air < 0):
397      !!----------------------------------------------------------------------
398
399            IF (s.lt.0. .and. ((stp.gt.-3.3).and.(stp.lt.0.))) THEN
400                fh = 0.1_wp+0.03_wp*stp+0.9_wp*exp(4.8_wp*stp)
401                fe = fh
402            ELSE IF (s.lt.0. .and. stp.le.-3.3) THEN
403                fh = 0._wp
404                fe = fh
405            ELSE                                       ! --- for unstable condition
406                fh = 1.0_wp+0.63_wp*sqrt(stp)
407                fe = fh
408            ENDIF
409
410      !!----------------------------------------------------------------------
411      !! --- calculate the coefficient CH,CE,CD
412      !!----------------------------------------------------------------------
413
414            IF (wndm(ji,jj) >= 0. .AND. wndm(ji,jj) <= 2.2)       THEN
415                kku=1
416            ELSE IF (wndm(ji,jj) > 2.2 .AND. wndm(ji,jj) <= 5.0)  THEN
417                kku=2
418            ELSE IF (wndm(ji,jj) > 5.0 .AND. wndm(ji,jj) <= 8.0)  THEN
419                kku=3
420            ELSE IF (wndm(ji,jj) > 8.0 .AND. wndm(ji,jj) <= 25.0) THEN
421                kku=4
422            ELSE IF (wndm(ji,jj) > 25.0 )                         THEN
423                kku=5
424            ENDIF
425
426            ch(ji,jj) = ( a_h(kku) + b_h(kku) * wndm(ji,jj) ** p_h(kku)      &
427                        + c_h(kku) * (wndm(ji,jj)-8 ) **2) * fh
428
429            ce(ji,jj) = ( a_e(kku) + b_e(kku) * wndm(ji,jj) ** p_e(kku)      &
430                        + c_e(kku) * (wndm(ji,jj)-8 ) **2) * fe
431
432            ch(ji,jj) = ch(ji,jj) / 1000.0
433            ce(ji,jj) = ce(ji,jj) / 1000.0
434
435            IF (wndm(ji,jj)<0.3) THEN
436               ch(ji,jj) = 1.3e-03 * fh
437               ce(ji,jj) = 1.5e-03 * fe
438            ELSE IF(wndm(ji,jj)>50.0) THEN
439               ch(ji,jj) = 1.25e-03 * fh
440               ce(ji,jj) = 1.30e-03 * fe
441            ENDIF
442
443      !!----------------------------------------------------------------------
444      !! --- calculates the SENSIBLE HEAT FLUX in MKS ( watt/m*m )
445      !!----------------------------------------------------------------------
446
447            HA(ji,jj) = rhom(ji,jj)*cp*ch(ji,jj)*wndm(ji,jj)*deltemp
448
449      !!----------------------------------------------------------------------
450      !! ------------------ (iv)  latent heat
451      !! --- calculates the LATENT HEAT FLUX  ( watt/m*m )
452      !! --- ELAT = L*rho*Ce*|V|*[qs(Ts)-qa(t2d)]
453      !!----------------------------------------------------------------------
454
455            esre  = shms(ji,jj) - shnow(ji,jj)   ! --- calculates the term : qs(Ta)-qa(t2d)
456
457            cseep = ce(ji,jj) * wndm(ji,jj) * esre     ! --- calculates the term : Ce*|V|*[qs(Ts)-qa(t2d)]
458
459            evap(ji,jj) = (cseep * rhom(ji,jj))  ! in [kg/m2/sec] !! --- calculates the EVAPORATION RATE [m/yr]
460
461            elat(ji,jj) = rhom(ji,jj) * cseep * heatlat(sst(ji,jj))
462
463      !!----------------------------------------------------------------------
464      !! --- calculates the Drag Coefficient
465      !!----------------------------------------------------------------------
466
467      !!----------------------------------------------------------------------
468      !! --- deltemp should be (Ts - Ta) in the formula estimating
469      !! --- drag coefficient
470      !!----------------------------------------------------------------------
471
472              IF( .NOT. ln_cdgw ) THEN
473                 cdx(ji,jj) = cd_HR(wndm(ji,jj),deltemp)
474              ENDIF
475
476          END DO
477      END DO
478
479      !!----------------------------------------------------------------------
480      !! --- calculates the wind stresses in MKS ( newton/m*m )
481      !! ---            taux= rho*Cd*|V|u      tauy= rho*Cd*|V|v
482      !!----------------------------------------------------------------------
483
484       taux(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windu(:,:)
485       tauy(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windv(:,:)
486
487      CALL wrk_dealloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
488      CALL wrk_dealloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
489      !
490      IF( nn_timing == 1 )  CALL timing_stop('fluxes_mfs')
491      !
492   END SUBROUTINE fluxes_mfs
493
494
495      REAL(wp) FUNCTION cd_HR(speed,delt)
496      !!----------------------------------------------------------------------
497      !! --- calculates the Drag Coefficient as a function of the abs. value
498      !! --- of the wind velocity ( Hellermann and Rosenstein )
499      !!----------------------------------------------------------------------
500
501       REAL(wp), INTENT(in) :: speed,delt
502       REAL(wp), PARAMETER  :: a1=0.934e-3 , a2=0.788e-4, a3=0.868e-4     
503       REAL(wp), PARAMETER  :: a4=-0.616e-6, a5=-.120e-5, a6=-.214e-5
504
505        cd_HR = a1 + a2*speed + a3*delt + a4*speed*speed        &
506           + a5*delt*delt  + a6*speed*delt
507
508      END FUNCTION cd_HR
509
510      REAL(wp) function HEATLAT(t)
511      !!----------------------------------------------------------------------
512      !! --- calculates the Latent Heat of Vaporization ( J/kg ) as function of
513      !! --- the temperature ( Celsius degrees )
514      !! --- ( from A. Gill  pag. 607 )
515      !!
516      !! --- Constant Latent Heat of Vaporization ( Rosati,Miyakoda 1988 )
517      !!     L = 2.501e+6  (MKS)
518      !!----------------------------------------------------------------------
519
520      REAL(wp) , intent(in) :: t
521
522      heatlat = 2.5008e+6 -2.3e+3*t
523
524      END FUNCTION HEATLAT
525
526
527   SUBROUTINE qshort(hour,alat,alon,cldnow,qsw)
528      !!----------------------------------------------------------------------
529      !!                    ***  ROUTINE qshort  ***
530      !!
531      !! ** Purpose :   Compute Solar Radiation
532      !!
533      !! ** Method  :   Compute Solar Radiation according Astronomical
534      !!                formulae
535      !!
536      !! References :   Reed RK (1975) and Reed RK (1977)
537      !!
538      !! Note: alat,alon - (lat, lon)  in radians
539      !!----------------------------------------------------------------------
540        REAL(wp), INTENT (in) :: hour
541
542        REAL(wp), INTENT(in ), DIMENSION(:,:) :: alat,alon
543        REAL(wp), INTENT(in ), DIMENSION(:,:) :: cldnow
544        REAL(wp), INTENT(out), DIMENSION(:,:) :: qsw
545        REAL(wp), DIMENSION(12) :: alpham
546
547        REAL(wp), PARAMETER ::   eclips=23.439* (3.141592653589793_wp / 180._wp)
548        REAL(wp), PARAMETER ::   solar = 1350.
549        REAL(wp), PARAMETER ::   tau = 0.7
550        REAL(wp), PARAMETER ::   aozone = 0.09
551        REAL(wp), PARAMETER ::   yrdays = 360.
552        REAL(wp) :: days, th0,th02,th03, sundec, thsun, coszen, qatten
553        REAL(wp) :: qzer, qdir,qdiff,qtot,tjul,sunbet
554        REAL(wp) :: albedo
555        INTEGER :: jj, ji
556
557      !!----------------------------------------------------------------------
558      !! --- albedo monthly values from Payne (1972) as means of the values
559      !! --- at 40N and 30N for the Atlantic Ocean ( hence the same latitudinal
560      !! --- band of the Mediterranean Sea ) :
561      !!----------------------------------------------------------------------
562
563        data alpham /0.095,0.08,0.065,0.065,0.06,0.06,0.06,0.06,        &
564                    0.065,0.075,0.09,0.10/
565
566      !!----------------------------------------------------------------------
567      !!   days is the number of days elapsed until the day=nday_year
568      !!----------------------------------------------------------------------
569        days = nday_year -1.
570        th0  = 2.*rpi*days/yrdays
571        th02 = 2.*th0
572        th03 = 3.*th0
573
574      !! --- sun declination :
575      !!----------------------------------------------------------------------
576        sundec = 0.006918 - 0.399912*cos(th0) + 0.070257*sin(th0) -   &
577                          0.006758*cos(th02) + 0.000907*sin(th02) -   &
578                          0.002697*cos(th03) + 0.001480*sin(th03)
579
580      DO jj = 1, jpj
581         DO ji = 1, jpi
582
583      !! --- sun hour angle :
584      !!----------------------------------------------------------------------
585          thsun = (hour -12.)*15.*rad + alon(ji,jj)
586
587      !! --- cosine of the solar zenith angle :
588      !!----------------------------------------------------------------------
589          coszen =sin(alat(ji,jj))*sin(sundec)                 &
590                    +cos(alat(ji,jj))*cos(sundec)*cos(thsun)
591
592          IF(coszen .LE. 5.035D-04) THEN
593            coszen = 0.0
594            qatten = 0.0
595          ELSE
596            qatten = tau**(1./coszen)
597          END IF
598
599          qzer  = coszen * solar *                                 &
600                  (1.0+1.67E-2*cos(rpi*2.*(days-3.0)/365.0))**2
601          qdir  = qzer * qatten
602          qdiff = ((1.-aozone)*qzer - qdir) * 0.5
603          qtot  =  qdir + qdiff
604          tjul = (days -81.)*rad
605
606      !! --- sin of the solar noon altitude in radians :
607      !!----------------------------------------------------------------------
608          sunbet=sin(alat(ji,jj))*sin(eclips*sin(tjul)) +   &
609                 cos(alat(ji,jj))*cos(eclips*sin(tjul))
610
611      !! --- solar noon altitude in degrees :
612      !!----------------------------------------------------------------------
613
614         sunbet = asin(sunbet)/rad
615
616      !!----------------------------------------------------------------------
617      !! --- calculates the albedo according to Payne (1972)
618      !!----------------------------------------------------------------------
619
620         albedo = alpham(nmonth)
621
622      !!----------------------------------------------------------------------
623      !! --- ( radiation as from Reed(1977), Simpson and Paulson(1979) )
624      !! --- calculates SHORT WAVE FLUX ( watt/m*m )
625      !! --- ( Rosati,Miyakoda 1988 ; eq. 3.8 )
626      !!----------------------------------------------------------------------
627
628          IF(cldnow(ji,jj).LT.0.3) THEN
629             qsw(ji,jj) = qtot * (1.-albedo)
630          ELSE
631             qsw(ji,jj) = qtot*(1.-0.62*cldnow(ji,jj)              &
632                                + .0019*sunbet)*(1.-albedo)
633          ENDIF
634
635         END DO
636      END DO
637
638   END SUBROUTINE qshort
639
640
641   !!======================================================================
642
643END MODULE sbcblk_mfs
Note: See TracBrowser for help on using the repository browser.