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 @ 5845

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

#1613: vvl by default: suppression of domzgr_substitute.h90

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