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_core.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

Last change on this file was 11101, checked in by frrh, 5 years ago

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