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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 6004

Last change on this file since 6004 was 6004, checked in by gm, 8 years ago

#1613: vvl by default, step III: Merge with the trunk (free surface simplification) (see wiki)
  • 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   USE sbcwave  ,ONLY :   cdn_wave !wave module
20   !
21   USE iom            ! I/O manager library
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! distribued memory computing library
24   USE wrk_nemo       ! work arrays
25   USE timing         ! Timing
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
27   USE prtctl         ! Print control
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   sbc_blk_mfs       ! routine called in sbcmod module
33     
34   INTEGER , PARAMETER ::   jpfld   = 7         ! maximum number of files to read
35   INTEGER , PARAMETER ::   jp_wndi = 1         ! index of 10m wind velocity (i-component) (m/s) at T-point
36   INTEGER , PARAMETER ::   jp_wndj = 2         ! index of 10m wind velocity (j-component) (m/s) at T-point
37   INTEGER , PARAMETER ::   jp_clc  = 3         ! index of total cloud cover               ( % )
38   INTEGER , PARAMETER ::   jp_msl  = 4         ! index of mean sea level pressure         (Pa)
39   INTEGER , PARAMETER ::   jp_tair = 5         ! index of 10m air temperature             (Kelvin)
40   INTEGER , PARAMETER ::   jp_rhm  = 6         ! index of dew point temperature           (Kelvin)
41   INTEGER , PARAMETER ::   jp_prec = 7         ! index of total precipitation (rain+snow) (Kg/m2/s)
42   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf ! structure of input fields (file informations, fields read)
43         
44   !! * Substitutions
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
48   !! $Id$
49   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE sbc_blk_mfs( kt )
54      !!---------------------------------------------------------------------
55      !!                    ***  ROUTINE sbc_blk_mfs  ***
56      !!                   
57      !! ** Purpose :   provide at each time step the surface ocean fluxes
58      !!      (momentum, heat, freshwater, runoff is added later in the code)
59      !!
60      !! ** Method  : (1) READ Atmospheric data from NetCDF files:
61      !!      the 10m wind velocity (i-component) (m/s) at T-point
62      !!      the 10m wind velocity (j-component) (m/s) at T-point
63      !!      the 2m Dew point Temperature        (k)
64      !!      the Cloud COver                     (%)
65      !!      the 2m air temperature              (Kelvin)
66      !!      the Mean Sea Level Preesure         (hPa)
67      !!      the Climatological Precipitation    (kg/m2/s)
68      !!              (2) CALL blk_oce_mfs
69      !!
70      !!      Computes:
71      !!      Solar Radiation using Reed formula (1975, 1977)
72      !!      Net Long wave radiation using Bignami et al. (1995)
73      !!      Latent and Sensible heat using Kondo (1975)
74      !!      Drag coeff using Hllerman and Rosenstein (1983)
75      !!      C A U T I O N : never mask the surface stress fields
76      !!                      the stress is assumed to be in the mesh referential
77      !!                      i.e. the (i,j) referential
78      !!
79      !! ** Action  :   defined at each time-step at the air-sea interface
80      !!              - utau, vtau  i- and j-component of the wind stress
81      !!              - taum        wind stress module at T-point
82      !!              - wndm        10m wind module at T-point over free ocean or leads in presence of sea-ice
83      !!              - qns, qsr    non-slor and solar heat flux
84      !!              - emp         evaporation minus precipitation
85      !!----------------------------------------------------------------------
86      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sh_now   ! specific humidity at T-point
87      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  catm     ! Cover
88      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  alonl    ! Lon
89      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  alatl    ! Lat
90      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gsst     ! SST
91      !!---------------------------------------------------------------------
92      !! Local fluxes variables
93      !!---------------------------------------------------------------------
94      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  qbw     ! Net Long wave
95      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ha      ! Sesnible
96      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  elat    ! Latent
97      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  evap    ! evaporation rate
98
99      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
100      !!
101      INTEGER  :: ierror                          ! return error code
102      INTEGER  :: ifpr     ! dummy loop indice
103      INTEGER  :: jj,ji    ! dummy loop arguments
104      INTEGER  ::   ios    ! Local integer output status for namelist read
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      IF( nn_timing == 1 )  CALL timing_start('sbc_blk_mfs')
124      !
125      !                                         ! ====================== !
126      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
127         !                                      ! ====================== !
128         ALLOCATE( sh_now(jpi,jpj), catm(jpi,jpj), alonl(jpi,jpj), alatl(jpi,jpj),     &
129         &        gsst(jpi,jpj),  qbw(jpi,jpj),    ha(jpi,jpj),  elat(jpi,jpj),     &
130         &        evap(jpi,jpj), STAT=ierror )
131
132         IF( ierror /= 0 )   CALL ctl_warn('sbc_blk_mfs: failed to allocate arrays')
133
134         REWIND( numnam_ref )              ! Namelist namsbc_msf in reference namelist : MFS files
135         READ  ( numnam_ref, namsbc_mfs, IOSTAT = ios, ERR = 901)
136901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_mfs in reference namelist', lwp )
137
138         REWIND( numnam_cfg )              ! Namelist namsbc_msf in configuration namelist : MFS files
139         READ  ( numnam_cfg, namsbc_mfs, IOSTAT = ios, ERR = 902 )
140902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_mfs in configuration namelist', lwp )
141         IF(lwm) WRITE ( numond, namsbc_mfs )
142         !
143         ! store namelist information in an array
144         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
145         slf_i(jp_clc ) = sn_clc    ;   slf_i(jp_msl ) = sn_msl
146         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_rhm)  = sn_rhm
147         slf_i(jp_prec) = sn_prec   ; 
148         !
149         ALLOCATE( sf(jpfld), STAT=ierror )         ! set sf structure
150         IF( ierror > 0 ) THEN
151            CALL ctl_stop( 'sbc_blk_mfs: unable to allocate sf structure' )   ;   RETURN
152         ENDIF
153         DO ifpr= 1, jpfld
154            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
155            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
156         END DO
157         ! fill sf with slf_i and control print
158         CALL fld_fill( sf, slf_i, cn_dir,'sbc_blk_mfs','bulk formulation for ocean SBC', 'namsbc_mfs' )
159            !
160      ENDIF
161
162         CALL fld_read( kt, nn_fsbc, sf )                   ! input fields provided at the current time-step
163
164         catm(:,:)   = 0.0    ! initializze cloud cover variable
165         sh_now(:,:) = 0.0    ! initializze specifif humidity variable
166
167         DO jj = 2, jpjm1
168            DO ji = fs_2, fs_jpim1   ! vector opt.
169
170         ! Calculate Specific Humidity
171         !-------------------------------------------------
172            sh_now(ji,jj) = (1/1.22) * onsea * par1 * EXP(par2/sf(jp_rhm)%fnow(ji,jj,1))
173
174         ! Normalize Clouds
175         !-------------------------------------------------
176            catm(ji,jj)   = max(0.0,min(1.0,sf(jp_clc)%fnow(ji,jj,1)*0.01))
177
178            END DO
179         END DO
180
181         ! wind module at 10m
182         !--------------------------------------------------
183         wndm(:,:) = SQRT(  sf(jp_wndi)%fnow(:,:,1) * sf(jp_wndi)%fnow(:,:,1)   &
184              &             + sf(jp_wndj)%fnow(:,:,1) * sf(jp_wndj)%fnow(:,:,1)  )
185
186         ! Some conv for fluxes computation
187         !-------------------------------------------------
188         alonl(:,:) = glamt(:,:) * rad
189         alatl(:,:) = gphit(:,:) * rad
190         gsst(:,:)  = tsn(:,:,1,jp_tem)  * tmask(:,:,1)
191
192         IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
193
194         ! Force to zero the output of fluxes
195         !-------------------------------------------------
196          qsr(:,:)  = 0.0 ; qbw(:,:)  = 0.0 ; 
197          ha(:,:)   = 0.0 ; elat(:,:) = 0.0 ; 
198          evap(:,:) = 0.0 ; utau(:,:) = 0.0 ; 
199          vtau(:,:) = 0.0
200
201          CALL lbc_lnk( sf(jp_wndi)%fnow(:,:,1), 'T', -1. )
202          CALL lbc_lnk( sf(jp_wndj)%fnow(:,:,1), 'T', -1. )
203
204          act_hour = (( nsec_year / rday ) - INT (nsec_year / rday)) * rjjhh
205
206          CALL fluxes_mfs(alatl,alonl,act_hour,                                &     ! input static
207                            gsst(:,:),sf(jp_tair)%fnow(:,:,1),sh_now(:,:),     &     ! input dynamic
208                            sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),  &     ! input dynamic
209                            sf(jp_msl)%fnow(:,:,1) , catm(:,:) ,               &     ! input dynamic
210                            qsr,qbw,ha,elat,evap,utau,vtau)                          ! output
211
212         ! Shortwave radiation
213         !--------------------------------------------------
214          qsr(:,:) = qsr(:,:) * tmask(:,:,1)
215
216         ! total non solar heat flux over water
217         !--------------------------------------------------
218          qns(:,:) = -1. * ( qbw(:,:) + ha(:,:) + elat(:,:) )
219          qns(:,:) = qns(:,:)*tmask(:,:,1)
220
221         ! mask the wind module at 10m
222         !--------------------------------------------------
223          wndm(:,:) = wndm(:,:) * tmask(:,:,1)
224
225         !   wind stress module (taum) into T-grid
226         !--------------------------------------------------
227          taum(:,:) = SQRT( utau(:,:) * utau(:,:) + vtau(:,:) * vtau(:,:) ) * tmask(:,:,1)
228
229          CALL lbc_lnk( taum, 'T', 1. )
230
231         ! Interpolate utau, vtau into the grid_V and grid_V
232         !-------------------------------------------------
233      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
234      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
235         DO jj = 1, jpjm1
236            DO ji = 1, fs_jpim1
237               utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( utau(ji,jj) * tmask(ji,jj,1) &
238               &                                + utau(ji+1,jj) * tmask(ji+1,jj,1) )        &
239               &                 * MAX(tmask(ji,jj,1),tmask(ji+1,jj  ,1))
240               vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vtau(ji,jj) * tmask(ji,jj,1) &
241               &                                + vtau(ji,jj+1) * tmask(ji,jj+1,1) )        &
242               &                 * MAX(tmask(ji,jj,1),tmask(ji  ,jj+1,1))
243            END DO
244         END DO
245
246         CALL lbc_lnk( utau(:,:), 'U', -1. )
247         CALL lbc_lnk( vtau(:,:), 'V', -1. )
248
249         ! for basin budget and cooerence
250         !--------------------------------------------------
251         emp(:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
252
253         CALL iom_put( "qlw_oce",   qbw  )                 ! output downward longwave heat over the ocean
254         CALL iom_put( "qsb_oce", - ha   )                 ! output downward sensible heat over the ocean
255         CALL iom_put( "qla_oce", - elat )                 ! output downward latent   heat over the ocean
256         CALL iom_put( "qns_oce",   qns  )                 ! output downward non solar heat over the ocean
257
258      ENDIF
259      !
260      IF( nn_timing == 1 )  CALL timing_stop('sbc_blk_mfs')
261      !
262   END SUBROUTINE sbc_blk_mfs
263 
264 
265   SUBROUTINE fluxes_mfs(alat,alon,hour,                               &
266        sst,tnow,shnow,unow,vnow,mslnow,cldnow,qsw,qbw,ha,elat,        &
267        evap,taux,tauy)
268      !!----------------------------------------------------------------------
269      !!                    ***  ROUTINE fluxes_mfs  ***
270      !!
271      !! --- it provides SURFACE HEAT and MOMENTUM FLUXES in MKS :
272      !!
273      !!  1)   Water flux (WFLUX)                 [ watt/m*m ]
274      !!  2)   Short wave flux (QSW)              [ watt/m*m ] Reed 1977
275      !!  3)   Long wave flux backward (QBW)      [ watt/m*m ]
276      !!  4)   Latent heat of evaporation (ELAT)  [ watt/m*m ]
277      !!  5)   Sensible heat flux   (HA)          [ watt/m*m ]
278      !!  6)   Wind stress x-component   (TAUX)   [ newton/m*m ]
279      !!  7)   Wind stress y-component   (TAUY)   [ newton/m*m ]
280      !!
281      !!----------------------------------------------------------------------
282      USE sbcblk_core, ONLY: turb_core_2z ! For wave coupling and Tair/rh from 2 to 10m
283
284      REAL(wp), INTENT(in   ) :: hour
285      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: sst, unow, alat , alon
286      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: vnow, cldnow, mslnow
287      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: qsw, qbw, ha, elat
288      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: evap,taux,tauy
289      REAL(wp), INTENT(inout), DIMENSION (:,:) :: tnow , shnow
290
291      INTEGER :: ji,jj 
292      REAL(wp)  :: wair, vtnow, ea, deltemp, s, stp , fh , fe
293      REAL(wp)  :: esre, cseep
294
295      REAL(wp), DIMENSION (:,:), POINTER ::   rspeed, sh10now, t10now, cdx, ce, shms
296      REAL(wp), DIMENSION (:,:), POINTER ::   rhom, sstk, ch, rel_windu, rel_windv
297      !!----------------------------------------------------------------------
298      !!     coefficients ( in MKS )  :
299      !!----------------------------------------------------------------------
300
301      REAL(wp), PARAMETER ::  ps = 1013.    ! --- surface air pressure
302      REAL(wp), PARAMETER ::  expsi=0.622   ! --- expsi
303      REAL(wp), PARAMETER ::  rd=287.       ! --- dry air gas constant
304      REAL(wp), PARAMETER ::  cp=1005.      ! --- specific heat capacity
305      REAL(wp), PARAMETER ::  onsea=0.98    ! --- specific humidity factors
306      REAL(wp), PARAMETER ::  par1=640380.  ! [Kg/m3]
307      REAL(wp), PARAMETER ::  par2=-5107.4  ! [K]
308
309      !---------------------------------------------------------------------
310      !--- define Kondo parameters
311      !---------------------------------------------------------------------
312
313      REAL(wp), DIMENSION(5) :: a_h = (/0.0,0.927,1.15,1.17,1.652/)
314      REAL(wp), DIMENSION(5) :: a_e = (/0.0,0.969,1.18,1.196,1.68/)
315      REAL(wp), DIMENSION(5) :: b_h = (/1.185,0.0546,0.01,0.0075,-0.017/)
316      REAL(wp), DIMENSION(5) :: b_e = (/1.23,0.0521,0.01,0.008,-0.016/)
317      REAL(wp), DIMENSION(5) :: c_h = (/0.0,0.0,0.0,-0.00045,0.0/)
318      REAL(wp), DIMENSION(5) :: c_e = (/0.0,0.0,0.0,-0.0004,0.0/)
319      REAL(wp), DIMENSION(5) :: p_h = (/-0.157,1.0,1.0,1.0,1.0/)
320      REAL(wp), DIMENSION(5) :: p_e = (/-0.16,1.0,1.0,1.0,1.0/)
321      INTEGER :: kku                        !index varing with wind speed
322      !
323      IF( nn_timing == 1 )  CALL timing_start('fluxes_mfs')
324      !
325      CALL wrk_alloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
326      CALL wrk_alloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
327
328      !!----------------------------------------------------------------------
329      !! ------------------ (i)      short wave
330      !!----------------------------------------------------------------------
331
332       CALL qshort(hour,alat,alon,cldnow,qsw)
333
334          rel_windu(:,:) = 0.0_wp
335          rel_windv(:,:) = 0.0_wp
336
337       DO jj = 2, jpj
338          DO ji = fs_2, jpi
339           rel_windu(ji,jj) = unow(ji,jj) - 0.5_wp * ( ssu_m(ji-1,jj) + ssu_m(ji,jj) )
340           rel_windv(ji,jj) = vnow(ji,jj) - 0.5_wp * ( ssv_m(ji,jj-1) + ssv_m(ji,jj) )
341          END DO
342       END DO
343
344       rspeed(:,:)= SQRT(rel_windu(:,:)*rel_windu(:,:)   &
345         &                   + rel_windv(:,:)*rel_windv(:,:)) 
346
347       sstk(:,:) = sst(:,:) + rtt                          !- SST data in Kelvin degrees
348       shms(:,:) = (1/1.22)*onsea*par1*EXP(par2/sstk(:,:)) !- Saturation Specific Humidity
349
350      ! --- Transport temperature and humidity from 2m to 10m
351      !----------------------------------------------------------------------
352
353      t10now(:,:) = 0.0           ;   sh10now(:,:)= 0.0
354      ! Note that air temp is converted in potential temp
355      CALL turb_core_2z(2.,10.,sstk,tnow+2*0.0098,shms,shnow,rspeed,        &
356         &              Cdx,Ch,Ce,t10now,sh10now )
357      tnow(:,:)  = t10now(:,:)    ;    shnow(:,:) = sh10now(:,:)
358
359      !!----------------------------------------------------------------------
360      !! ------------------ (ii)    net long wave
361      !!----------------------------------------------------------------------
362
363      DO jj = 1, jpj
364         DO ji = 1, jpi
365            wair = shnow(ji,jj) / (1 - shnow(ji,jj))    ! mixing ratio of the air (Wallace and Hobbs)
366            vtnow = (tnow(ji,jj)*(expsi+wair))/(expsi*(1.+wair))   ! virtual temperature of air
367            rhom(ji,jj) = 100.*(ps/rd)/vtnow                       ! density of the moist air
368            ea   = (wair  / (wair  + 0.622 )) * mslnow(ji,jj)
369
370            qbw(ji,jj) = emic*stefan*( sstk(ji,jj)**4. )                    &
371                 - ( stefan*( tnow(ji,jj)**4. ) * ( 0.653 + 0.00535*ea ) )  &
372                   * ( 1. + 0.1762*( cldnow(ji,jj)**2. ) )
373
374         END DO
375      END DO
376
377      DO jj = 1, jpj
378         DO ji = 1, jpi
379      !!----------------------------------------------------------------------
380      !! ------------------ (iii)   sensible heat
381      !!----------------------------------------------------------------------
382
383      !! --- calculates the term :      ( Ts - Ta )
384      !!----------------------------------------------------------------------
385            deltemp = sstk(ji,jj) - tnow (ji,jj)
386
387      !!----------------------------------------------------------------------
388      !! --- variable turbulent exchange coefficients ( from Kondo 1975 )
389      !! --- calculate the Neutral Transfer Coefficent using an empiric formula
390      !! --- by Kondo et al. Then it applies the diabatic approximation.
391      !!----------------------------------------------------------------------
392
393            s = deltemp/(wndm(ji,jj)**2.)   !! --- calculate S
394            stp = s*abs(s)/(abs(s)+0.01)    !! --- calculate the Stability Parameter
395
396      !!----------------------------------------------------------------------
397      !! --- for stable condition (sst-t_air < 0):
398      !!----------------------------------------------------------------------
399
400            IF (s.lt.0. .and. ((stp.gt.-3.3).and.(stp.lt.0.))) THEN
401                fh = 0.1_wp+0.03_wp*stp+0.9_wp*exp(4.8_wp*stp)
402                fe = fh
403            ELSE IF (s.lt.0. .and. stp.le.-3.3) THEN
404                fh = 0._wp
405                fe = fh
406            ELSE                                       ! --- for unstable condition
407                fh = 1.0_wp+0.63_wp*sqrt(stp)
408                fe = fh
409            ENDIF
410
411      !!----------------------------------------------------------------------
412      !! --- calculate the coefficient CH,CE,CD
413      !!----------------------------------------------------------------------
414
415            IF (wndm(ji,jj) >= 0. .AND. wndm(ji,jj) <= 2.2)       THEN
416                kku=1
417            ELSE IF (wndm(ji,jj) > 2.2 .AND. wndm(ji,jj) <= 5.0)  THEN
418                kku=2
419            ELSE IF (wndm(ji,jj) > 5.0 .AND. wndm(ji,jj) <= 8.0)  THEN
420                kku=3
421            ELSE IF (wndm(ji,jj) > 8.0 .AND. wndm(ji,jj) <= 25.0) THEN
422                kku=4
423            ELSE IF (wndm(ji,jj) > 25.0 )                         THEN
424                kku=5
425            ENDIF
426
427            ch(ji,jj) = ( a_h(kku) + b_h(kku) * wndm(ji,jj) ** p_h(kku)      &
428                        + c_h(kku) * (wndm(ji,jj)-8 ) **2) * fh
429
430            ce(ji,jj) = ( a_e(kku) + b_e(kku) * wndm(ji,jj) ** p_e(kku)      &
431                        + c_e(kku) * (wndm(ji,jj)-8 ) **2) * fe
432
433            ch(ji,jj) = ch(ji,jj) / 1000.0
434            ce(ji,jj) = ce(ji,jj) / 1000.0
435
436            IF (wndm(ji,jj)<0.3) THEN
437               ch(ji,jj) = 1.3e-03 * fh
438               ce(ji,jj) = 1.5e-03 * fe
439            ELSE IF(wndm(ji,jj)>50.0) THEN
440               ch(ji,jj) = 1.25e-03 * fh
441               ce(ji,jj) = 1.30e-03 * fe
442            ENDIF
443
444      !!----------------------------------------------------------------------
445      !! --- calculates the SENSIBLE HEAT FLUX in MKS ( watt/m*m )
446      !!----------------------------------------------------------------------
447
448            HA(ji,jj) = rhom(ji,jj)*cp*ch(ji,jj)*wndm(ji,jj)*deltemp
449
450      !!----------------------------------------------------------------------
451      !! ------------------ (iv)  latent heat
452      !! --- calculates the LATENT HEAT FLUX  ( watt/m*m )
453      !! --- ELAT = L*rho*Ce*|V|*[qs(Ts)-qa(t2d)]
454      !!----------------------------------------------------------------------
455
456            esre  = shms(ji,jj) - shnow(ji,jj)   ! --- calculates the term : qs(Ta)-qa(t2d)
457
458            cseep = ce(ji,jj) * wndm(ji,jj) * esre     ! --- calculates the term : Ce*|V|*[qs(Ts)-qa(t2d)]
459
460            evap(ji,jj) = (cseep * rhom(ji,jj))  ! in [kg/m2/sec] !! --- calculates the EVAPORATION RATE [m/yr]
461
462            elat(ji,jj) = rhom(ji,jj) * cseep * heatlat(sst(ji,jj))
463
464      !!----------------------------------------------------------------------
465      !! --- calculates the Drag Coefficient
466      !!----------------------------------------------------------------------
467
468      !!----------------------------------------------------------------------
469      !! --- deltemp should be (Ts - Ta) in the formula estimating
470      !! --- drag coefficient
471      !!----------------------------------------------------------------------
472
473              IF( .NOT. ln_cdgw ) THEN
474                 cdx(ji,jj) = cd_HR(wndm(ji,jj),deltemp)
475              ENDIF
476
477          END DO
478      END DO
479
480      !!----------------------------------------------------------------------
481      !! --- calculates the wind stresses in MKS ( newton/m*m )
482      !! ---            taux= rho*Cd*|V|u      tauy= rho*Cd*|V|v
483      !!----------------------------------------------------------------------
484
485       taux(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windu(:,:)
486       tauy(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windv(:,:)
487
488      CALL wrk_dealloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
489      CALL wrk_dealloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
490      !
491      IF( nn_timing == 1 )  CALL timing_stop('fluxes_mfs')
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.