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

source: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 3107

Last change on this file since 3107 was 3107, checked in by poddo, 12 years ago

correct bug to to tsn in mfs bulk

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