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

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 8145

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

#1909 bug fix to allow reproducibility when sbcblk_mfs is used

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