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

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90

Last change on this file was 11277, checked in by kingr, 5 years ago

Merged Juan's changes for running AMM15 woth wave coupling.
Corrected minor logic error to allow AMM7-uncoupled to reproduce earlier results.
Few line spacing changes to allow merging with OBS br and trunk rvn 5518.

File size: 29.9 KB
RevLine 
[3105]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
[3294]22   USE wrk_nemo        ! work arrays
23   USE timing          ! Timing
[3105]24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE prtctl          ! Print control
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   sbc_blk_mfs       ! routine called in sbcmod module
31     
32   INTEGER , PARAMETER ::   jpfld   = 7         ! maximum number of files to read
33   INTEGER , PARAMETER ::   jp_wndi = 1         ! index of 10m wind velocity (i-component) (m/s) at T-point
34   INTEGER , PARAMETER ::   jp_wndj = 2         ! index of 10m wind velocity (j-component) (m/s) at T-point
35   INTEGER , PARAMETER ::   jp_clc  = 3         ! index of total cloud cover               ( % )
36   INTEGER , PARAMETER ::   jp_msl  = 4         ! index of mean sea level pressure         (Pa)
37   INTEGER , PARAMETER ::   jp_tair = 5         ! index of 10m air temperature             (Kelvin)
38   INTEGER , PARAMETER ::   jp_rhm  = 6         ! index of dew point temperature           (Kelvin)
39   INTEGER , PARAMETER ::   jp_prec = 7         ! index of total precipitation (rain+snow) (Kg/m2/s)
40   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf ! structure of input fields (file informations, fields read)
41         
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
[5215]47   !! $Id$
[3105]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
[4990]83      !!              - wndm        10m wind module at T-point over free ocean or leads in presence of sea-ice
[3105]84      !!              - qns, qsr    non-slor and solar heat flux
[3625]85      !!              - emp         evaporation minus precipitation
[3105]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
[4147]105      INTEGER  ::   ios    ! Local integer output status for namelist read
[3105]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      !!---------------------------------------------------------------------
[3294]123      !
124      IF( nn_timing == 1 )  CALL timing_start('sbc_blk_mfs')
125      !
[3105]126      !                                         ! ====================== !
127      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
128         !                                      ! ====================== !
[4147]129         ALLOCATE( sh_now(jpi,jpj), catm(jpi,jpj), alonl(jpi,jpj), alatl(jpi,jpj),     &
[3105]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')
[4147]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 )
[4624]142         IF(lwm) WRITE ( numond, namsbc_mfs )
[4147]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' )
[3105]160            !
[4147]161      ENDIF
[3105]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
[3107]191         gsst(:,:)  = tsn(:,:,1,jp_tem)  * tmask(:,:,1)
[3105]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         !-------------------------------------------------
[4990]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
[3105]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) &
[4990]239               &                                + utau(ji+1,jj) * tmask(ji+1,jj,1) )        &
240               &                 * MAX(tmask(ji,jj,1),tmask(ji+1,jj  ,1))
[3105]241               vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vtau(ji,jj) * tmask(ji,jj,1) &
[4990]242               &                                + vtau(ji,jj+1) * tmask(ji,jj+1,1) )        &
243               &                 * MAX(tmask(ji,jj,1),tmask(ji  ,jj+1,1))
[3105]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!CDIR COLLAPSE
253           emp (:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
254!CDIR COLLAPSE
255
256         CALL iom_put( "qlw_oce",   qbw  )                 ! output downward longwave heat over the ocean
257         CALL iom_put( "qsb_oce", - ha   )                 ! output downward sensible heat over the ocean
258         CALL iom_put( "qla_oce", - elat )                 ! output downward latent   heat over the ocean
259         CALL iom_put( "qns_oce",   qns  )                 ! output downward non solar heat over the ocean
260
261      ENDIF
[3294]262      !
263      IF( nn_timing == 1 )  CALL timing_stop('sbc_blk_mfs')
264      !
[3105]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      USE sbcblk_core, ONLY: turb_core_2z ! For wave coupling and Tair/rh from 2 to 10m
286
287      REAL(wp), INTENT(in   ) :: hour
288      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: sst, unow, alat , alon
289      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: vnow, cldnow, mslnow
290      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: qsw, qbw, ha, elat
291      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: evap,taux,tauy
292      REAL(wp), INTENT(inout), DIMENSION (:,:) :: tnow , shnow
293
294      INTEGER :: ji,jj 
295      REAL(wp)  :: wair, vtnow, ea, deltemp, s, stp , fh , fe
296      REAL(wp)  :: esre, cseep
297
[3294]298      REAL(wp), DIMENSION (:,:), POINTER ::   rspeed, sh10now, t10now, cdx, ce, shms
299      REAL(wp), DIMENSION (:,:), POINTER ::   rhom, sstk, ch, rel_windu, rel_windv
[3105]300      !!----------------------------------------------------------------------
301      !!     coefficients ( in MKS )  :
302      !!----------------------------------------------------------------------
303
304      REAL(wp), PARAMETER ::  ps = 1013.    ! --- surface air pressure
305      REAL(wp), PARAMETER ::  expsi=0.622   ! --- expsi
306      REAL(wp), PARAMETER ::  rd=287.       ! --- dry air gas constant
307      REAL(wp), PARAMETER ::  cp=1005.      ! --- specific heat capacity
308      REAL(wp), PARAMETER ::  onsea=0.98    ! --- specific humidity factors
309      REAL(wp), PARAMETER ::  par1=640380.  ! [Kg/m3]
310      REAL(wp), PARAMETER ::  par2=-5107.4  ! [K]
311
312      !---------------------------------------------------------------------
313      !--- define Kondo parameters
314      !---------------------------------------------------------------------
315
316      REAL(wp), DIMENSION(5) :: a_h = (/0.0,0.927,1.15,1.17,1.652/)
317      REAL(wp), DIMENSION(5) :: a_e = (/0.0,0.969,1.18,1.196,1.68/)
318      REAL(wp), DIMENSION(5) :: b_h = (/1.185,0.0546,0.01,0.0075,-0.017/)
319      REAL(wp), DIMENSION(5) :: b_e = (/1.23,0.0521,0.01,0.008,-0.016/)
320      REAL(wp), DIMENSION(5) :: c_h = (/0.0,0.0,0.0,-0.00045,0.0/)
321      REAL(wp), DIMENSION(5) :: c_e = (/0.0,0.0,0.0,-0.0004,0.0/)
322      REAL(wp), DIMENSION(5) :: p_h = (/-0.157,1.0,1.0,1.0,1.0/)
323      REAL(wp), DIMENSION(5) :: p_e = (/-0.16,1.0,1.0,1.0,1.0/)
324      INTEGER :: kku                        !index varing with wind speed
[3294]325      !
326      IF( nn_timing == 1 )  CALL timing_start('fluxes_mfs')
327      !
328      CALL wrk_alloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
329      CALL wrk_alloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
[3105]330
331      !!----------------------------------------------------------------------
332      !! ------------------ (i)      short wave
333      !!----------------------------------------------------------------------
334
335       CALL qshort(hour,alat,alon,cldnow,qsw)
336
337          rel_windu(:,:) = 0.0_wp
338          rel_windv(:,:) = 0.0_wp
339
340       DO jj = 2, jpj
341          DO ji = fs_2, jpi
342           rel_windu(ji,jj) = unow(ji,jj) - 0.5_wp * ( ssu_m(ji-1,jj) + ssu_m(ji,jj) )
343           rel_windv(ji,jj) = vnow(ji,jj) - 0.5_wp * ( ssv_m(ji,jj-1) + ssv_m(ji,jj) )
344          END DO
345       END DO
346
347       rspeed(:,:)= SQRT(rel_windu(:,:)*rel_windu(:,:)   &
348         &                   + rel_windv(:,:)*rel_windv(:,:)) 
349
350       sstk(:,:) = sst(:,:) + rtt                          !- SST data in Kelvin degrees
351       shms(:,:) = (1/1.22)*onsea*par1*EXP(par2/sstk(:,:)) !- Saturation Specific Humidity
352
353      ! --- Transport temperature and humidity from 2m to 10m
354      !----------------------------------------------------------------------
355
356      t10now(:,:) = 0.0           ;   sh10now(:,:)= 0.0
357      ! Note that air temp is converted in potential temp
358      CALL turb_core_2z(2.,10.,sstk,tnow+2*0.0098,shms,shnow,rspeed,        &
359         &              Cdx,Ch,Ce,t10now,sh10now )
360      tnow(:,:)  = t10now(:,:)    ;    shnow(:,:) = sh10now(:,:)
361
362      !!----------------------------------------------------------------------
363      !! ------------------ (ii)    net long wave
364      !!----------------------------------------------------------------------
365
366      DO jj = 1, jpj
367         DO ji = 1, jpi
368            wair = shnow(ji,jj) / (1 - shnow(ji,jj))    ! mixing ratio of the air (Wallace and Hobbs)
369            vtnow = (tnow(ji,jj)*(expsi+wair))/(expsi*(1.+wair))   ! virtual temperature of air
370            rhom(ji,jj) = 100.*(ps/rd)/vtnow                       ! density of the moist air
371            ea   = (wair  / (wair  + 0.622 )) * mslnow(ji,jj)
372
373            qbw(ji,jj) = emic*stefan*( sstk(ji,jj)**4. )                    &
374                 - ( stefan*( tnow(ji,jj)**4. ) * ( 0.653 + 0.00535*ea ) )  &
375                   * ( 1. + 0.1762*( cldnow(ji,jj)**2. ) )
376
377         END DO
378      END DO
379
380      DO jj = 1, jpj
381         DO ji = 1, jpi
382      !!----------------------------------------------------------------------
383      !! ------------------ (iii)   sensible heat
384      !!----------------------------------------------------------------------
385
386      !! --- calculates the term :      ( Ts - Ta )
387      !!----------------------------------------------------------------------
388            deltemp = sstk(ji,jj) - tnow (ji,jj)
389
390      !!----------------------------------------------------------------------
391      !! --- variable turbulent exchange coefficients ( from Kondo 1975 )
392      !! --- calculate the Neutral Transfer Coefficent using an empiric formula
393      !! --- by Kondo et al. Then it applies the diabatic approximation.
394      !!----------------------------------------------------------------------
395
396            s = deltemp/(wndm(ji,jj)**2.)   !! --- calculate S
397            stp = s*abs(s)/(abs(s)+0.01)    !! --- calculate the Stability Parameter
398
399      !!----------------------------------------------------------------------
400      !! --- for stable condition (sst-t_air < 0):
401      !!----------------------------------------------------------------------
402
403            IF (s.lt.0. .and. ((stp.gt.-3.3).and.(stp.lt.0.))) THEN
404                fh = 0.1_wp+0.03_wp*stp+0.9_wp*exp(4.8_wp*stp)
405                fe = fh
406            ELSE IF (s.lt.0. .and. stp.le.-3.3) THEN
407                fh = 0._wp
408                fe = fh
409            ELSE                                       ! --- for unstable condition
410                fh = 1.0_wp+0.63_wp*sqrt(stp)
411                fe = fh
412            ENDIF
413
414      !!----------------------------------------------------------------------
415      !! --- calculate the coefficient CH,CE,CD
416      !!----------------------------------------------------------------------
417
418            IF (wndm(ji,jj) >= 0. .AND. wndm(ji,jj) <= 2.2)       THEN
419                kku=1
420            ELSE IF (wndm(ji,jj) > 2.2 .AND. wndm(ji,jj) <= 5.0)  THEN
421                kku=2
422            ELSE IF (wndm(ji,jj) > 5.0 .AND. wndm(ji,jj) <= 8.0)  THEN
423                kku=3
424            ELSE IF (wndm(ji,jj) > 8.0 .AND. wndm(ji,jj) <= 25.0) THEN
425                kku=4
426            ELSE IF (wndm(ji,jj) > 25.0 )                         THEN
427                kku=5
428            ENDIF
429
430            ch(ji,jj) = ( a_h(kku) + b_h(kku) * wndm(ji,jj) ** p_h(kku)      &
431                        + c_h(kku) * (wndm(ji,jj)-8 ) **2) * fh
432
433            ce(ji,jj) = ( a_e(kku) + b_e(kku) * wndm(ji,jj) ** p_e(kku)      &
434                        + c_e(kku) * (wndm(ji,jj)-8 ) **2) * fe
435
436            ch(ji,jj) = ch(ji,jj) / 1000.0
437            ce(ji,jj) = ce(ji,jj) / 1000.0
438
439            IF (wndm(ji,jj)<0.3) THEN
440               ch(ji,jj) = 1.3e-03 * fh
441               ce(ji,jj) = 1.5e-03 * fe
442            ELSE IF(wndm(ji,jj)>50.0) THEN
443               ch(ji,jj) = 1.25e-03 * fh
444               ce(ji,jj) = 1.30e-03 * fe
445            ENDIF
446
447      !!----------------------------------------------------------------------
448      !! --- calculates the SENSIBLE HEAT FLUX in MKS ( watt/m*m )
449      !!----------------------------------------------------------------------
450
451            HA(ji,jj) = rhom(ji,jj)*cp*ch(ji,jj)*wndm(ji,jj)*deltemp
452
453      !!----------------------------------------------------------------------
454      !! ------------------ (iv)  latent heat
455      !! --- calculates the LATENT HEAT FLUX  ( watt/m*m )
456      !! --- ELAT = L*rho*Ce*|V|*[qs(Ts)-qa(t2d)]
457      !!----------------------------------------------------------------------
458
459            esre  = shms(ji,jj) - shnow(ji,jj)   ! --- calculates the term : qs(Ta)-qa(t2d)
460
461            cseep = ce(ji,jj) * wndm(ji,jj) * esre     ! --- calculates the term : Ce*|V|*[qs(Ts)-qa(t2d)]
462
463            evap(ji,jj) = (cseep * rhom(ji,jj))  ! in [kg/m2/sec] !! --- calculates the EVAPORATION RATE [m/yr]
464
465            elat(ji,jj) = rhom(ji,jj) * cseep * heatlat(sst(ji,jj))
466
467      !!----------------------------------------------------------------------
468      !! --- calculates the Drag Coefficient
469      !!----------------------------------------------------------------------
470
471      !!----------------------------------------------------------------------
472      !! --- deltemp should be (Ts - Ta) in the formula estimating
473      !! --- drag coefficient
474      !!----------------------------------------------------------------------
475
476              IF( .NOT. ln_cdgw ) THEN
477                 cdx(ji,jj) = cd_HR(wndm(ji,jj),deltemp)
478              ENDIF
479
480          END DO
481      END DO
482
483      !!----------------------------------------------------------------------
484      !! --- calculates the wind stresses in MKS ( newton/m*m )
485      !! ---            taux= rho*Cd*|V|u      tauy= rho*Cd*|V|v
486      !!----------------------------------------------------------------------
487
488       taux(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windu(:,:)
489       tauy(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windv(:,:)
490
[3294]491      CALL wrk_dealloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
492      CALL wrk_dealloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
493      !
494      IF( nn_timing == 1 )  CALL timing_stop('fluxes_mfs')
495      !
[3105]496   END SUBROUTINE fluxes_mfs
497
498
499      REAL(wp) FUNCTION cd_HR(speed,delt)
500      !!----------------------------------------------------------------------
501      !! --- calculates the Drag Coefficient as a function of the abs. value
502      !! --- of the wind velocity ( Hellermann and Rosenstein )
503      !!----------------------------------------------------------------------
504
505       REAL(wp), INTENT(in) :: speed,delt
506       REAL(wp), PARAMETER  :: a1=0.934e-3 , a2=0.788e-4, a3=0.868e-4     
507       REAL(wp), PARAMETER  :: a4=-0.616e-6, a5=-.120e-5, a6=-.214e-5
508
509        cd_HR = a1 + a2*speed + a3*delt + a4*speed*speed        &
510           + a5*delt*delt  + a6*speed*delt
511
512      END FUNCTION cd_HR
513
514      REAL(wp) function HEATLAT(t)
515      !!----------------------------------------------------------------------
516      !! --- calculates the Latent Heat of Vaporization ( J/kg ) as function of
517      !! --- the temperature ( Celsius degrees )
518      !! --- ( from A. Gill  pag. 607 )
519      !!
520      !! --- Constant Latent Heat of Vaporization ( Rosati,Miyakoda 1988 )
521      !!     L = 2.501e+6  (MKS)
522      !!----------------------------------------------------------------------
523
524      REAL(wp) , intent(in) :: t
525
526      heatlat = 2.5008e+6 -2.3e+3*t
527
528      END FUNCTION HEATLAT
529
530
531   SUBROUTINE qshort(hour,alat,alon,cldnow,qsw)
532      !!----------------------------------------------------------------------
533      !!                    ***  ROUTINE qshort  ***
534      !!
535      !! ** Purpose :   Compute Solar Radiation
536      !!
537      !! ** Method  :   Compute Solar Radiation according Astronomical
538      !!                formulae
539      !!
540      !! References :   Reed RK (1975) and Reed RK (1977)
541      !!
542      !! Note: alat,alon - (lat, lon)  in radians
543      !!----------------------------------------------------------------------
544        REAL(wp), INTENT (in) :: hour
545
546        REAL(wp), INTENT(in ), DIMENSION(:,:) :: alat,alon
547        REAL(wp), INTENT(in ), DIMENSION(:,:) :: cldnow
548        REAL(wp), INTENT(out), DIMENSION(:,:) :: qsw
549        REAL(wp), DIMENSION(12) :: alpham
550
551        REAL(wp), PARAMETER ::   eclips=23.439* (3.141592653589793_wp / 180._wp)
552        REAL(wp), PARAMETER ::   solar = 1350.
553        REAL(wp), PARAMETER ::   tau = 0.7
554        REAL(wp), PARAMETER ::   aozone = 0.09
555        REAL(wp), PARAMETER ::   yrdays = 360.
556        REAL(wp) :: days, th0,th02,th03, sundec, thsun, coszen, qatten
557        REAL(wp) :: qzer, qdir,qdiff,qtot,tjul,sunbet
558        REAL(wp) :: albedo
559        INTEGER :: jj, ji
560
561      !!----------------------------------------------------------------------
562      !! --- albedo monthly values from Payne (1972) as means of the values
563      !! --- at 40N and 30N for the Atlantic Ocean ( hence the same latitudinal
564      !! --- band of the Mediterranean Sea ) :
565      !!----------------------------------------------------------------------
566
567        data alpham /0.095,0.08,0.065,0.065,0.06,0.06,0.06,0.06,        &
568                    0.065,0.075,0.09,0.10/
569
570      !!----------------------------------------------------------------------
571      !!   days is the number of days elapsed until the day=nday_year
572      !!----------------------------------------------------------------------
573        days = nday_year -1.
574        th0  = 2.*rpi*days/yrdays
575        th02 = 2.*th0
576        th03 = 3.*th0
577
578      !! --- sun declination :
579      !!----------------------------------------------------------------------
580        sundec = 0.006918 - 0.399912*cos(th0) + 0.070257*sin(th0) -   &
581                          0.006758*cos(th02) + 0.000907*sin(th02) -   &
582                          0.002697*cos(th03) + 0.001480*sin(th03)
583
584      DO jj = 1, jpj
585         DO ji = 1, jpi
586
587      !! --- sun hour angle :
588      !!----------------------------------------------------------------------
589          thsun = (hour -12.)*15.*rad + alon(ji,jj)
590
591      !! --- cosine of the solar zenith angle :
592      !!----------------------------------------------------------------------
593          coszen =sin(alat(ji,jj))*sin(sundec)                 &
594                    +cos(alat(ji,jj))*cos(sundec)*cos(thsun)
595
596          IF(coszen .LE. 5.035D-04) THEN
597            coszen = 0.0
598            qatten = 0.0
599          ELSE
600            qatten = tau**(1./coszen)
601          END IF
602
603          qzer  = coszen * solar *                                 &
604                  (1.0+1.67E-2*cos(rpi*2.*(days-3.0)/365.0))**2
605          qdir  = qzer * qatten
606          qdiff = ((1.-aozone)*qzer - qdir) * 0.5
607          qtot  =  qdir + qdiff
608          tjul = (days -81.)*rad
609
610      !! --- sin of the solar noon altitude in radians :
611      !!----------------------------------------------------------------------
612          sunbet=sin(alat(ji,jj))*sin(eclips*sin(tjul)) +   &
613                 cos(alat(ji,jj))*cos(eclips*sin(tjul))
614
615      !! --- solar noon altitude in degrees :
616      !!----------------------------------------------------------------------
617
618         sunbet = asin(sunbet)/rad
619
620      !!----------------------------------------------------------------------
621      !! --- calculates the albedo according to Payne (1972)
622      !!----------------------------------------------------------------------
623
624         albedo = alpham(nmonth)
625
626      !!----------------------------------------------------------------------
627      !! --- ( radiation as from Reed(1977), Simpson and Paulson(1979) )
628      !! --- calculates SHORT WAVE FLUX ( watt/m*m )
629      !! --- ( Rosati,Miyakoda 1988 ; eq. 3.8 )
630      !!----------------------------------------------------------------------
631
632          IF(cldnow(ji,jj).LT.0.3) THEN
633             qsw(ji,jj) = qtot * (1.-albedo)
634          ELSE
635             qsw(ji,jj) = qtot*(1.-0.62*cldnow(ji,jj)              &
636                                + .0019*sunbet)*(1.-albedo)
637          ENDIF
638
639         END DO
640      END DO
641
642   END SUBROUTINE qshort
643
644
645   !!======================================================================
646
647END MODULE sbcblk_mfs
Note: See TracBrowser for help on using the repository browser.