source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 @ 11442

Last change on this file since 11442 was 11442, checked in by mattmartin, 18 months ago

Introduction of stochastic physics in NEMO, based on Andrea Storto's code.
For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.

File size: 51.5 KB
Line 
1MODULE sbcblk_core
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk_core  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!=====================================================================
6   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original code
7   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier) additions:
8   !!                           -  new bulk routine for efficiency
9   !!                           -  WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files !!!!
10   !!                           -  file names and file characteristics in namelist
11   !!                           -  Implement reading of 6-hourly fields
12   !!            3.0  !  2006-06  (G. Madec) sbc rewritting
13   !!             -   !  2006-12  (L. Brodeau) Original code for turb_core_2z
14   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
15   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
16   !!            3.4  !  2011-11  (C. Harris) Fill arrays required by CICE
17   !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   sbc_blk_core    : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea)
22   !!   blk_oce_core    : computes momentum, heat and freshwater fluxes over ocean
23   !!   blk_ice_core    : computes momentum, heat and freshwater fluxes over ice
24   !!   turb_core_2z    : Computes turbulent transfert coefficients
25   !!   cd_neutral_10m  : Estimate of the neutral drag coefficient at 10m
26   !!   psi_m           : universal profile stability function for momentum
27   !!   psi_h           : universal profile stability function for temperature and humidity
28   !!----------------------------------------------------------------------
29   USE oce             ! ocean dynamics and tracers
30   USE dom_oce         ! ocean space and time domain
31   USE phycst          ! physical constants
32   USE fldread         ! read input fields
33   USE sbc_oce         ! Surface boundary condition: ocean fields
34   USE cyclone         ! Cyclone 10m wind form trac of cyclone centres
35   USE sbcdcy          ! surface boundary condition: diurnal cycle
36   USE iom             ! I/O manager library
37   USE in_out_manager  ! I/O manager
38   USE lib_mpp         ! distribued memory computing library
39   USE wrk_nemo        ! work arrays
40   USE timing          ! Timing
41   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
42   USE prtctl          ! Print control
43   USE sbcwave, ONLY   :  cdn_wave ! wave module
44   USE sbc_ice         ! Surface boundary condition: ice fields
45   USE lib_fortran     ! to use key_nosignedzero
46#if defined key_lim3
47   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b
48   USE limthd_dh       ! for CALL lim_thd_snwblow
49#elif defined key_lim2
50   USE ice_2, ONLY     : u_ice, v_ice
51   USE par_ice_2
52#endif
53   USE stopack
54
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC   sbc_blk_core         ! routine called in sbcmod module
59#if defined key_lim2 || defined key_lim3
60   PUBLIC   blk_ice_core_tau     ! routine called in sbc_ice_lim module
61   PUBLIC   blk_ice_core_flx     ! routine called in sbc_ice_lim module
62#endif
63   PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module
64
65   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read
66   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
67   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
68   INTEGER , PARAMETER ::   jp_humi = 3           ! index of specific humidity               ( % )
69   INTEGER , PARAMETER ::   jp_qsr  = 4           ! index of solar heat                      (W/m2)
70   INTEGER , PARAMETER ::   jp_qlw  = 5           ! index of Long wave                       (W/m2)
71   INTEGER , PARAMETER ::   jp_tair = 6           ! index of 10m air temperature             (Kelvin)
72   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
73   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
74   INTEGER , PARAMETER ::   jp_tdif = 9           ! index of tau diff associated to HF tau   (N/m2)   at T-point
75
76   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
77
78   !                                             !!! CORE bulk parameters
79   REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density
80   REAL(wp), PARAMETER ::   cpa  = 1000.5         ! specific heat of air
81   REAL(wp), PARAMETER ::   Lv   =    2.5e6       ! latent heat of vaporization
82   REAL(wp), PARAMETER ::   Ls   =    2.839e6     ! latent heat of sublimation
83   REAL(wp), PARAMETER ::   Stef =    5.67e-8     ! Stefan Boltzmann constant
84   REAL(wp), PARAMETER ::   Cice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice
85   REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant
86
87   !                                  !!* Namelist namsbc_core : CORE bulk parameters
88   LOGICAL  ::   ln_taudif   ! logical flag to use the "mean of stress module - module of mean stress" data
89   REAL(wp) ::   rn_pfac     ! multiplication factor for precipitation
90   REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem)
91   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
92   REAL(wp), ALLOCATABLE, SAVE ::   rn_vfac0(:,:) ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
93   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
94   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
95   REAL(wp), PUBLIC :: rn_sfac ! multiplication factor for snow precipitation over sea-ice
96
97   !! * Substitutions
98#  include "domzgr_substitute.h90"
99#  include "vectopt_loop_substitute.h90"
100   !!----------------------------------------------------------------------
101   !! NEMO/OPA 3.7 , NEMO-consortium (2014)
102   !! $Id$
103   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
104   !!----------------------------------------------------------------------
105CONTAINS
106
107   SUBROUTINE sbc_blk_core( kt )
108      !!---------------------------------------------------------------------
109      !!                    ***  ROUTINE sbc_blk_core  ***
110      !!
111      !! ** Purpose :   provide at each time step the surface ocean fluxes
112      !!      (momentum, heat, freshwater and runoff)
113      !!
114      !! ** Method  : (1) READ each fluxes in NetCDF files:
115      !!      the 10m wind velocity (i-component) (m/s)    at T-point
116      !!      the 10m wind velocity (j-component) (m/s)    at T-point
117      !!      the 10m or 2m specific humidity     ( % )
118      !!      the solar heat                      (W/m2)
119      !!      the Long wave                       (W/m2)
120      !!      the 10m or 2m air temperature       (Kelvin)
121      !!      the total precipitation (rain+snow) (Kg/m2/s)
122      !!      the snow (solid prcipitation)       (kg/m2/s)
123      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T)
124      !!              (2) CALL blk_oce_core
125      !!
126      !!      C A U T I O N : never mask the surface stress fields
127      !!                      the stress is assumed to be in the (i,j) mesh referential
128      !!
129      !! ** Action  :   defined at each time-step at the air-sea interface
130      !!              - utau, vtau  i- and j-component of the wind stress
131      !!              - taum        wind stress module at T-point
132      !!              - wndm        wind speed  module at T-point over free ocean or leads in presence of sea-ice
133      !!              - qns, qsr    non-solar and solar heat fluxes
134      !!              - emp         upward mass flux (evapo. - precip.)
135      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present)
136      !!                            (set in limsbc(_2).F90)
137      !!
138      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
139      !!                   Brodeau et al. Ocean Modelling 2010
140      !!----------------------------------------------------------------------
141      INTEGER, INTENT(in) ::   kt   ! ocean time step
142      !
143      INTEGER  ::   ierror   ! return error code
144      INTEGER  ::   ifpr     ! dummy loop indice
145      INTEGER  ::   jfld     ! dummy loop arguments
146      INTEGER  ::   ios      ! Local integer output status for namelist read
147      !
148      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files
149      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
150      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
151      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 "
152      TYPE(FLD_N) ::   sn_tdif                                 !   "                                 "
153      NAMELIST/namsbc_core/ cn_dir , ln_taudif, rn_pfac, rn_efac, rn_vfac,  &
154         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
155         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
156         &                  sn_tdif, rn_zqt,  rn_zu, rn_sfac
157      !!---------------------------------------------------------------------
158      !
159      !                                         ! ====================== !
160      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
161         !                                      ! ====================== !
162         !
163         rn_sfac = 1._wp       ! Default to one if missing from namelist
164         REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters
165         READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901)
166901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in reference namelist', lwp )
167         !
168         REWIND( numnam_cfg )              ! Namelist namsbc_core in configuration namelist : CORE bulk parameters
169         READ  ( numnam_cfg, namsbc_core, IOSTAT = ios, ERR = 902 )
170902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in configuration namelist', lwp )
171
172         IF(lwm) WRITE( numond, namsbc_core )
173         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
174         IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   &
175            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' )
176         IF( ln_dm2dc .AND. sn_qsr%ln_tint ) THEN
177            CALL ctl_warn( 'sbc_blk_core: ln_dm2dc is taking care of the temporal interpolation of daily qsr',   &
178               &         '              ==> We force time interpolation = .false. for qsr' )
179            sn_qsr%ln_tint = .false.
180         ENDIF
181         !                                         ! store namelist information in an array
182         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
183         slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
184         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
185         slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
186         slf_i(jp_tdif) = sn_tdif
187         !
188         lhftau = ln_taudif                        ! do we use HF tau information?
189         jfld = jpfld - COUNT( (/.NOT. lhftau/) )
190         !
191         ALLOCATE( sf(jfld), STAT=ierror )         ! set sf structure
192         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core: unable to allocate sf structure' )
193         DO ifpr= 1, jfld
194            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
195            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
196         END DO
197         !                                         ! fill sf with slf_i and control print
198         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
199         !
200         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
201         !
202         ALLOCATE ( rn_vfac0(jpi,jpj) )
203         rn_vfac0(:,:) = rn_vfac
204         !
205      ENDIF
206
207      IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN
208         rn_vfac0(:,:) = rn_vfac
209         CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw )
210      ENDIF
211
212      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
213
214      !                                            ! compute the surface ocean fluxes using CORE bulk formulea
215      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( kt, sf, sst_m, ssu_m, ssv_m )
216
217#if defined key_cice
218      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
219         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
220         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
221         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
222         ENDIF
223         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)         
224         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
225         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
226         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
227         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
228         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
229      ENDIF
230#endif
231      !
232   END SUBROUTINE sbc_blk_core
233   
234   
235   SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv )
236      !!---------------------------------------------------------------------
237      !!                     ***  ROUTINE blk_core  ***
238      !!
239      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
240      !!      the ocean surface at each time step
241      !!
242      !! ** Method  :   CORE bulk formulea for the ocean using atmospheric
243      !!      fields read in sbc_read
244      !!
245      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
246      !!              - vtau    : j-component of the stress at V-point  (N/m2)
247      !!              - taum    : Wind stress module at T-point         (N/m2)
248      !!              - wndm    : Wind speed module at T-point          (m/s)
249      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
250      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
251      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
252      !!
253      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
254      !!---------------------------------------------------------------------
255      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
256      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
257      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
258      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
259      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
260      !
261      INTEGER  ::   ji, jj               ! dummy loop indices
262      REAL(wp) ::   zcoef_qsatw, zztmp   ! local variable
263      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point
264      REAL(wp), DIMENSION(:,:), POINTER ::   zqsatw            ! specific humidity at pst
265      REAL(wp), DIMENSION(:,:), POINTER ::   zqlw, zqsb        ! long wave and sensible heat fluxes
266      REAL(wp), DIMENSION(:,:), POINTER ::   zqla, zevap       ! latent heat fluxes and evaporation
267      REAL(wp), DIMENSION(:,:), POINTER ::   Cd                ! transfer coefficient for momentum      (tau)
268      REAL(wp), DIMENSION(:,:), POINTER ::   Ch                ! transfer coefficient for sensible heat (Q_sens)
269      REAL(wp), DIMENSION(:,:), POINTER ::   Ce                ! tansfert coefficient for evaporation   (Q_lat)
270      REAL(wp), DIMENSION(:,:), POINTER ::   zst               ! surface temperature in Kelvin
271      REAL(wp), DIMENSION(:,:), POINTER ::   zt_zu             ! air temperature at wind speed height
272      REAL(wp), DIMENSION(:,:), POINTER ::   zq_zu             ! air spec. hum.  at wind speed height
273      !!---------------------------------------------------------------------
274      !
275      IF( nn_timing == 1 )  CALL timing_start('blk_oce_core')
276      !
277      CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap )
278      CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
279      !
280      ! local scalars ( place there for vector optimisation purposes)
281      zcoef_qsatw = 0.98 * 640380. / rhoa
282     
283      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
284
285      ! ----------------------------------------------------------------------------- !
286      !      0   Wind components and module at T-point relative to the moving ocean   !
287      ! ----------------------------------------------------------------------------- !
288
289      ! ... components ( U10m - U_oce ) at T-point (unmasked)
290      zwnd_i(:,:) = 0.e0 
291      zwnd_j(:,:) = 0.e0
292#if defined key_cyclone
293      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
294      DO jj = 2, jpjm1
295         DO ji = fs_2, fs_jpim1   ! vect. opt.
296            sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
297            sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
298         END DO
299      END DO
300#endif
301      DO jj = 2, jpjm1
302         DO ji = fs_2, fs_jpim1   ! vect. opt.
303            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
304            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
305         END DO
306      END DO
307      CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
308      CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
309      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
310      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
311         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
312
313      ! ----------------------------------------------------------------------------- !
314      !      I   Radiative FLUXES                                                     !
315      ! ----------------------------------------------------------------------------- !
316
317      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
318      zztmp = 1. - albo
319      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
320      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
321      ENDIF
322
323      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
324      ! ----------------------------------------------------------------------------- !
325      !     II    Turbulent FLUXES                                                    !
326      ! ----------------------------------------------------------------------------- !
327
328      ! ... specific humidity at SST and IST
329      zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) )
330
331      ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point :
332      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   &
333         &               Cd, Ch, Ce, zt_zu, zq_zu )
334   
335      ! ... tau module, i and j component
336      DO jj = 1, jpj
337         DO ji = 1, jpi
338            zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj)
339            taum  (ji,jj) = zztmp * wndm  (ji,jj)
340            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
341            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
342         END DO
343      END DO
344
345      ! ... add the HF tau contribution to the wind stress module?
346      IF( lhftau ) THEN
347         taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
348      ENDIF
349      CALL iom_put( "taum_oce", taum )   ! output wind stress module
350
351      ! ... utau, vtau at U- and V_points, resp.
352      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
353      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
354      DO jj = 1, jpjm1
355         DO ji = 1, fs_jpim1
356            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
357               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
358            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
359               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
360         END DO
361      END DO
362      CALL lbc_lnk( utau(:,:), 'U', -1. )
363      CALL lbc_lnk( vtau(:,:), 'V', -1. )
364
365   
366      !  Turbulent fluxes over ocean
367      ! -----------------------------
368      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
369         !! q_air and t_air are (or "are almost") given at 10m (wind reference height)
370         zevap(:,:) = rn_efac*MAX( 0._wp,     rhoa*Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) )*wndm(:,:) ) ! Evaporation
371         zqsb (:,:) =                     cpa*rhoa*Ch(:,:)*( zst   (:,:) - sf(jp_tair)%fnow(:,:,1) )*wndm(:,:)   ! Sensible Heat
372      ELSE
373         !! q_air and t_air are not given at 10m (wind reference height)
374         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
375         zevap(:,:) = rn_efac*MAX( 0._wp,     rhoa*Ce(:,:)*( zqsatw(:,:) - zq_zu(:,:) )*wndm(:,:) )   ! Evaporation
376         zqsb (:,:) =                     cpa*rhoa*Ch(:,:)*( zst   (:,:) - zt_zu(:,:) )*wndm(:,:)     ! Sensible Heat
377      ENDIF
378      zqla (:,:) = Lv * zevap(:,:)                                                              ! Latent Heat
379
380      IF(ln_ctl) THEN
381         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce_core: zqla   : ', tab2d_2=Ce , clinfo2=' Ce  : ' )
382         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce_core: zqsb   : ', tab2d_2=Ch , clinfo2=' Ch  : ' )
383         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce_core: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
384         CALL prt_ctl( tab2d_1=zqsatw, clinfo1=' blk_oce_core: zqsatw : ', tab2d_2=zst, clinfo2=' zst : ' )
385         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce_core: utau   : ', mask1=umask,   &
386            &          tab2d_2=vtau  , clinfo2=              ' vtau : '  , mask2=vmask )
387         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_core: wndm   : ')
388         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce_core: zst    : ')
389      ENDIF
390       
391      ! ----------------------------------------------------------------------------- !
392      !     III    Total FLUXES                                                       !
393      ! ----------------------------------------------------------------------------- !
394      !
395      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
396         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
397      !
398      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
399         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
400         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
401         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
402         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
403         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
404         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1)
405      !
406#if defined key_lim3
407      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by LIM3)
408      qsr_oce(:,:) = qsr(:,:)
409#endif
410      !
411      IF ( nn_ice == 0 ) THEN
412         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
413         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
414         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
415         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
416         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
417         CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
418         CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
419         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac   ! output total precipitation [kg/m2/s]
420         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac   ! output solid precipitation [kg/m2/s]
421         CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow
422         CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation
423      ENDIF
424      !
425      IF(ln_ctl) THEN
426         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce_core: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
427         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_core: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
428         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce_core: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
429         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce_core: utau   : ', mask1=umask,   &
430            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
431      ENDIF
432      !
433      CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap )
434      CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
435      !
436      IF( nn_timing == 1 )  CALL timing_stop('blk_oce_core')
437      !
438   END SUBROUTINE blk_oce_core
439 
440   
441#if defined key_lim2 || defined key_lim3
442   SUBROUTINE blk_ice_core_tau
443      !!---------------------------------------------------------------------
444      !!                     ***  ROUTINE blk_ice_core_tau  ***
445      !!
446      !! ** Purpose :   provide the surface boundary condition over sea-ice
447      !!
448      !! ** Method  :   compute momentum using CORE bulk
449      !!                formulea, ice variables and read atmospheric fields.
450      !!                NB: ice drag coefficient is assumed to be a constant
451      !!---------------------------------------------------------------------
452      INTEGER  ::   ji, jj    ! dummy loop indices
453      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2
454      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point
455      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point
456      !!---------------------------------------------------------------------
457      !
458      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau')
459      !
460      ! local scalars ( place there for vector optimisation purposes)
461      zcoef_wnorm  = rhoa * Cice
462      zcoef_wnorm2 = rhoa * Cice * 0.5
463
464!!gm brutal....
465      utau_ice  (:,:) = 0._wp
466      vtau_ice  (:,:) = 0._wp
467      wndm_ice  (:,:) = 0._wp
468!!gm end
469
470      ! ----------------------------------------------------------------------------- !
471      !    Wind components and module relative to the moving ocean ( U10m - U_ice )   !
472      ! ----------------------------------------------------------------------------- !
473      SELECT CASE( cp_ice_msh )
474      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
475         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked)
476         DO jj = 2, jpjm1
477            DO ji = 2, jpim1   ! B grid : NO vector opt
478               ! ... scalar wind at I-point (fld being at T-point)
479               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)    &
480                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) &
481                  &      - rn_vfac0(ji,jj) * u_ice(ji,jj)
482               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)    &
483                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) &
484                  &      - rn_vfac0(ji,jj) * v_ice(ji,jj)
485               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
486               ! ... ice stress at I-point
487               utau_ice(ji,jj) = zwnorm_f * zwndi_f
488               vtau_ice(ji,jj) = zwnorm_f * zwndj_f
489               ! ... scalar wind at T-point (fld being at T-point)
490               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1)                                         &
491                  &      - rn_vfac0(ji,jj) * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   &
492                  &                                  + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  )
493               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1)                                         &
494                  &      - rn_vfac0(ji,jj) * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   &
495                  &                                  + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  )
496               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
497            END DO
498         END DO
499         CALL lbc_lnk( utau_ice, 'I', -1. )
500         CALL lbc_lnk( vtau_ice, 'I', -1. )
501         CALL lbc_lnk( wndm_ice, 'T',  1. )
502         !
503      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean)
504         DO jj = 2, jpj
505            DO ji = fs_2, jpi   ! vect. opt.
506               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
507               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
508               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
509            END DO
510         END DO
511         DO jj = 2, jpjm1
512            DO ji = fs_2, fs_jpim1   ! vect. opt.
513               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )             &
514                  &              * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) &
515                  &                  - rn_vfac0(ji,jj) * u_ice(ji,jj) )
516               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )             &
517                  &              * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) &
518                  &                  - rn_vfac0(ji,jj) * v_ice(ji,jj) )
519            END DO
520         END DO
521         CALL lbc_lnk( utau_ice, 'U', -1. )
522         CALL lbc_lnk( vtau_ice, 'V', -1. )
523         CALL lbc_lnk( wndm_ice, 'T',  1. )
524         !
525      END SELECT
526
527      IF(ln_ctl) THEN
528         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice_core: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
529         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice_core: wndm_ice : ')
530      ENDIF
531
532      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau')
533     
534   END SUBROUTINE blk_ice_core_tau
535
536
537   SUBROUTINE blk_ice_core_flx( ptsu, palb )
538      !!---------------------------------------------------------------------
539      !!                     ***  ROUTINE blk_ice_core_flx  ***
540      !!
541      !! ** Purpose :   provide the surface boundary condition over sea-ice
542      !!
543      !! ** Method  :   compute heat and freshwater exchanged
544      !!                between atmosphere and sea-ice using CORE bulk
545      !!                formulea, ice variables and read atmmospheric fields.
546      !!
547      !! caution : the net upward water flux has with mm/day unit
548      !!---------------------------------------------------------------------
549      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu          ! sea ice surface temperature
550      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb          ! ice albedo (all skies)
551      !!
552      INTEGER  ::   ji, jj, jl    ! dummy loop indices
553      REAL(wp) ::   zst2, zst3
554      REAL(wp) ::   zcoef_dqlw, zcoef_dqla, zcoef_dqsb
555      REAL(wp) ::   zztmp, z1_lsub                               ! temporary variable
556      !!
557      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice
558      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice
559      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice
560      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice
561      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw       ! evaporation and snw distribution after wind blowing (LIM3)
562      !!---------------------------------------------------------------------
563      !
564      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx')
565      !
566      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
567
568      ! local scalars ( place there for vector optimisation purposes)
569      zcoef_dqlw   = 4.0 * 0.95 * Stef
570      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8)
571      zcoef_dqsb   = rhoa * cpa * Cice
572
573      zztmp = 1. / ( 1. - albo )
574      !                                     ! ========================== !
575      DO jl = 1, jpl                        !  Loop over ice categories  !
576         !                                  ! ========================== !
577         DO jj = 1 , jpj
578            DO ji = 1, jpi
579               ! ----------------------------!
580               !      I   Radiative FLUXES   !
581               ! ----------------------------!
582               zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
583               zst3 = ptsu(ji,jj,jl) * zst2
584               ! Short Wave (sw)
585               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
586               ! Long  Wave (lw)
587               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
588               ! lw sensitivity
589               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                               
590
591               ! ----------------------------!
592               !     II    Turbulent FLUXES  !
593               ! ----------------------------!
594
595               ! ... turbulent heat fluxes
596               ! Sensible Heat
597               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
598               ! Latent Heat
599               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                           
600                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
601              ! Latent heat sensitivity for ice (Dqla/Dt)
602               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
603                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) )
604               ELSE
605                  dqla_ice(ji,jj,jl) = 0._wp
606               ENDIF
607
608               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
609               z_dqsb(ji,jj,jl) = zcoef_dqsb * wndm_ice(ji,jj)
610
611               ! ----------------------------!
612               !     III    Total FLUXES     !
613               ! ----------------------------!
614               ! Downward Non Solar flux
615               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
616               ! Total non solar heat flux sensitivity for ice
617               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
618            END DO
619            !
620         END DO
621         !
622      END DO
623      !
624      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s]
625      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s]
626      CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation
627      CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation
628
629#if defined  key_lim3
630      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
631
632      ! --- evaporation --- !
633      z1_lsub = 1._wp / Lsub
634      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation
635      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT
636      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean
637
638      ! --- evaporation minus precipitation --- !
639      zsnw(:,:) = 0._wp
640      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing
641      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
642      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
643      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
644
645      ! --- heat flux associated with emp --- !
646      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst
647         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
648         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
649         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
650      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
651         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
652
653      ! --- total solar and non solar fluxes --- !
654      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
655      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
656
657      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
658      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
659
660      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
661      DO jl = 1, jpl
662         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
663                                   ! But we do not have Tice => consider it at 0 degC => evap=0
664      END DO
665
666      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
667#endif
668
669      !--------------------------------------------------------------------
670      ! FRACTIONs of net shortwave radiation which is not absorbed in the
671      ! thin surface layer and penetrates inside the ice cover
672      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
673      !
674      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
675      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
676      !
677      !
678      IF(ln_ctl) THEN
679         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice_core: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
680         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
681         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
682         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice_core: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
683         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice_core: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
684         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice_core: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
685      ENDIF
686
687      CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )
688      !
689      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx')
690     
691   END SUBROUTINE blk_ice_core_flx
692#endif
693
694   SUBROUTINE turb_core_2z( zt, zu, sst, T_zt, q_sat, q_zt, dU,    &
695      &                      Cd, Ch, Ce , T_zu, q_zu )
696      !!----------------------------------------------------------------------
697      !!                      ***  ROUTINE  turb_core  ***
698      !!
699      !! ** Purpose :   Computes turbulent transfert coefficients of surface
700      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
701      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
702      !!
703      !! ** Method : Monin Obukhov Similarity Theory
704      !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)
705      !!
706      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
707      !!
708      !! ** Last update: Laurent Brodeau, June 2014:
709      !!    - handles both cases zt=zu and zt/=zu
710      !!    - optimized: less 2D arrays allocated and less operations
711      !!    - better first guess of stability by checking air-sea difference of virtual temperature
712      !!       rather than temperature difference only...
713      !!    - added function "cd_neutral_10m" that uses the improved parametrization of
714      !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!
715      !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them
716      !!      => 'vkarmn' and 'grav'
717      !!----------------------------------------------------------------------
718      REAL(wp), INTENT(in   )                     ::   zt       ! height for T_zt and q_zt                   [m]
719      REAL(wp), INTENT(in   )                     ::   zu       ! height for dU                              [m]
720      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature              [Kelvin]
721      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   T_zt     ! potential air temperature            [Kelvin]
722      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_sat    ! sea surface specific humidity         [kg/kg]
723      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                 [kg/kg]
724      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   dU       ! relative wind module at zu            [m/s]
725      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
726      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
727      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
728      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   T_zu     ! air temp. shifted at zu                     [K]
729      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. hum.  shifted at zu               [kg/kg]
730      !
731      INTEGER ::   j_itt
732      INTEGER , PARAMETER ::   nb_itt = 5       ! number of itterations
733      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at different height than U
734      !
735      REAL(wp), DIMENSION(:,:), POINTER ::   U_zu          ! relative wind at zu                            [m/s]
736      REAL(wp), DIMENSION(:,:), POINTER ::   Ce_n10        ! 10m neutral latent coefficient
737      REAL(wp), DIMENSION(:,:), POINTER ::   Ch_n10        ! 10m neutral sensible coefficient
738      REAL(wp), DIMENSION(:,:), POINTER ::   sqrt_Cd_n10   ! root square of Cd_n10
739      REAL(wp), DIMENSION(:,:), POINTER ::   sqrt_Cd       ! root square of Cd
740      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_u        ! stability parameter at height zu
741      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_t        ! stability parameter at height zt
742      REAL(wp), DIMENSION(:,:), POINTER ::   zpsi_h_u, zpsi_m_u
743      REAL(wp), DIMENSION(:,:), POINTER ::   ztmp0, ztmp1, ztmp2
744      REAL(wp), DIMENSION(:,:), POINTER ::   stab          ! 1st stability test integer
745      !!----------------------------------------------------------------------
746
747      IF( nn_timing == 1 )  CALL timing_start('turb_core_2z')
748   
749      CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd )
750      CALL wrk_alloc( jpi,jpj, zeta_u, stab )
751      CALL wrk_alloc( jpi,jpj, zpsi_h_u, zpsi_m_u, ztmp0, ztmp1, ztmp2 )
752
753      l_zt_equal_zu = .FALSE.
754      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
755
756      IF( .NOT. l_zt_equal_zu )   CALL wrk_alloc( jpi,jpj, zeta_t )
757
758      U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
759
760      !! First guess of stability:
761      ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt
762      stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
763
764      !! Neutral coefficients at 10m:
765      IF( ln_cdgw ) THEN      ! wave drag case
766         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
767         ztmp0   (:,:) = cdn_wave(:,:)
768      ELSE
769         ztmp0 = cd_neutral_10m( U_zu )
770      ENDIF
771      sqrt_Cd_n10 = SQRT( ztmp0 )
772      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 )
773      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))
774   
775      !! Initializing transf. coeff. with their first guess neutral equivalents :
776      Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10
777
778      !! Initializing values at z_u with z_t values:
779      T_zu = T_zt   ;   q_zu = q_zt
780
781      !!  * Now starting iteration loop
782      DO j_itt=1, nb_itt
783         !
784         ztmp1 = T_zu - sst   ! Updating air/sea differences
785         ztmp2 = q_zu - q_sat 
786
787         ! Updating turbulent scales :   (L&Y 2004 eq. (7))
788         ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta*
789         ztmp2  = Ce/sqrt_Cd*ztmp2    ! q*
790       
791         ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu
792
793         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
794         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 
795         !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu)
796
797         !! Stability parameters :
798         zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u )
799         zpsi_h_u = psi_h( zeta_u )
800         zpsi_m_u = psi_m( zeta_u )
801       
802         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
803         IF ( .NOT. l_zt_equal_zu ) THEN
804            zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),10.0), zeta_t )
805            stab = LOG(zu/zt) - zpsi_h_u + psi_h(zeta_t)  ! stab just used as temp array!!!
806            T_zu = T_zt + ztmp1/vkarmn*stab    ! ztmp1 is still theta*
807            q_zu = q_zt + ztmp2/vkarmn*stab    ! ztmp2 is still q*
808            q_zu = max(0., q_zu)
809         END IF
810       
811         IF( ln_cdgw ) THEN      ! surface wave case
812            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
813            Cd      = sqrt_Cd * sqrt_Cd
814         ELSE
815           ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
816           !   In very rare low-wind conditions, the old way of estimating the
817           !   neutral wind speed at 10m leads to a negative value that causes the code
818           !   to crash. To prevent this a threshold of 0.25m/s is imposed.
819           ztmp0 = MAX( 0.25 , U_zu/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)) ) !  U_n10
820           ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10
821           sqrt_Cd_n10 = sqrt(ztmp0)
822       
823           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b)
824           stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability
825           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)
826
827           !! Update of transfer coefficients:
828           ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a)
829           Cd      = ztmp0 / ( ztmp1*ztmp1 )   
830           sqrt_Cd = SQRT( Cd )
831         ENDIF
832         !
833         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
834         ztmp2 = sqrt_Cd / sqrt_Cd_n10
835         ztmp1 = 1. + Ch_n10*ztmp0               
836         Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
837         !
838         ztmp1 = 1. + Ce_n10*ztmp0               
839         Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
840         !
841      END DO
842
843      CALL wrk_dealloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd )
844      CALL wrk_dealloc( jpi,jpj, zeta_u, stab )
845      CALL wrk_dealloc( jpi,jpj, zpsi_h_u, zpsi_m_u, ztmp0, ztmp1, ztmp2 )
846
847      IF( .NOT. l_zt_equal_zu ) CALL wrk_dealloc( jpi,jpj, zeta_t )
848
849      IF( nn_timing == 1 )  CALL timing_stop('turb_core_2z')
850      !
851   END SUBROUTINE turb_core_2z
852
853
854   FUNCTION cd_neutral_10m( zw10 )
855      !!----------------------------------------------------------------------
856      !! Estimate of the neutral drag coefficient at 10m as a function
857      !! of neutral wind  speed at 10m
858      !!
859      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)
860      !!
861      !! Author: L. Brodeau, june 2014
862      !!----------------------------------------------------------------------   
863      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s)
864      REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m
865      !
866      REAL(wp), DIMENSION(:,:), POINTER ::   rgt33
867      !!----------------------------------------------------------------------   
868      !
869      CALL wrk_alloc( jpi,jpj, rgt33 )
870      !
871      !! When wind speed > 33 m/s => Cyclone conditions => special treatment
872      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1 
873      cd_neutral_10m = 1.e-3 * ( &
874         &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33.
875         &      + rgt33         *      2.34   )                                                    ! zw10 >= 33.
876      !
877      CALL wrk_dealloc( jpi,jpj, rgt33)
878      !
879   END FUNCTION cd_neutral_10m
880
881
882   FUNCTION psi_m(pta)   !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
883      !-------------------------------------------------------------------------------
884      ! universal profile stability function for momentum
885      !-------------------------------------------------------------------------------
886      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta
887      !
888      REAL(wp), DIMENSION(jpi,jpj)             :: psi_m
889      REAL(wp), DIMENSION(:,:), POINTER        :: X2, X, stabit
890      !-------------------------------------------------------------------------------
891      !
892      CALL wrk_alloc( jpi,jpj, X2, X, stabit )
893      !
894      X2 = SQRT( ABS( 1. - 16.*pta ) )  ;  X2 = MAX( X2 , 1. )   ;   X = SQRT( X2 )
895      stabit = 0.5 + SIGN( 0.5 , pta )
896      psi_m = -5.*pta*stabit  &                                                          ! Stable
897         &    + (1. - stabit)*(2.*LOG((1. + X)*0.5) + LOG((1. + X2)*0.5) - 2.*ATAN(X) + rpi*0.5)  ! Unstable
898      !
899      CALL wrk_dealloc( jpi,jpj, X2, X, stabit )
900      !
901   END FUNCTION psi_m
902
903
904   FUNCTION psi_h( pta )    !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
905      !-------------------------------------------------------------------------------
906      ! universal profile stability function for temperature and humidity
907      !-------------------------------------------------------------------------------
908      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pta
909      !
910      REAL(wp), DIMENSION(jpi,jpj)             ::   psi_h
911      REAL(wp), DIMENSION(:,:), POINTER        ::   X2, X, stabit
912      !-------------------------------------------------------------------------------
913      !
914      CALL wrk_alloc( jpi,jpj, X2, X, stabit )
915      !
916      X2 = SQRT( ABS( 1. - 16.*pta ) )   ;   X2 = MAX( X2 , 1. )   ;   X = SQRT( X2 )
917      stabit = 0.5 + SIGN( 0.5 , pta )
918      psi_h = -5.*pta*stabit   &                                       ! Stable
919         &    + (1. - stabit)*(2.*LOG( (1. + X2)*0.5 ))                ! Unstable
920      !
921      CALL wrk_dealloc( jpi,jpj, X2, X, stabit )
922      !
923   END FUNCTION psi_h
924
925   !!======================================================================
926END MODULE sbcblk_core
Note: See TracBrowser for help on using the repository browser.