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/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 3556

Last change on this file since 3556 was 3294, checked in by rblod, 12 years ago

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