source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 10 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo

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