source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 7 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.

  • Property svn:keywords set to Id
File size: 58.4 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   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   sbc_blk_core  : bulk formulation as ocean surface boundary condition
21   !!                   (forced mode, CORE bulk formulea)
22   !!   blk_oce_core  : ocean: computes momentum, heat and freshwater fluxes
23   !!   blk_ice_core  : ice  : computes momentum, heat and freshwater fluxes
24   !!   turb_core     : computes the CORE turbulent transfer coefficients
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE phycst          ! physical constants
29   USE fldread         ! read input fields
30   USE sbc_oce         ! Surface boundary condition: ocean fields
31   USE cyclone         ! Cyclone 10m wind form trac of cyclone centres
32   USE sbcdcy          ! surface boundary condition: diurnal cycle
33   USE iom             ! I/O manager library
34   USE in_out_manager  ! I/O manager
35   USE lib_mpp         ! distribued memory computing library
36   USE wrk_nemo        ! work arrays
37   USE timing          ! Timing
38   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
39   USE prtctl          ! Print control
40   USE sbcwave,ONLY :  cdn_wave !wave module
41#if defined key_lim3 || defined key_cice
42   USE sbc_ice         ! Surface boundary condition: ice fields
43#endif
44   USE lib_fortran     ! to use key_nosignedzero
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC   sbc_blk_core         ! routine called in sbcmod module
50   PUBLIC   blk_ice_core         ! routine called in sbc_ice_lim module
51   PUBLIC   blk_ice_meanqsr      ! routine called in sbc_ice_lim module
52   PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module
53
54   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read
55   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
56   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
57   INTEGER , PARAMETER ::   jp_humi = 3           ! index of specific humidity               ( % )
58   INTEGER , PARAMETER ::   jp_qsr  = 4           ! index of solar heat                      (W/m2)
59   INTEGER , PARAMETER ::   jp_qlw  = 5           ! index of Long wave                       (W/m2)
60   INTEGER , PARAMETER ::   jp_tair = 6           ! index of 10m air temperature             (Kelvin)
61   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
62   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
63   INTEGER , PARAMETER ::   jp_tdif = 9           ! index of tau diff associated to HF tau   (N/m2)   at T-point
64   
65   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
66         
67   !                                             !!! CORE bulk parameters
68   REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density
69   REAL(wp), PARAMETER ::   cpa  = 1000.5         ! specific heat of air
70   REAL(wp), PARAMETER ::   Lv   =    2.5e6       ! latent heat of vaporization
71   REAL(wp), PARAMETER ::   Ls   =    2.839e6     ! latent heat of sublimation
72   REAL(wp), PARAMETER ::   Stef =    5.67e-8     ! Stefan Boltzmann constant
73   REAL(wp), PARAMETER ::   Cice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice
74   REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant
75
76   !                                  !!* Namelist namsbc_core : CORE bulk parameters
77   LOGICAL  ::   ln_2m       ! logical flag for height of air temp. and hum
78   LOGICAL  ::   ln_taudif   ! logical flag to use the "mean of stress module - module of mean stress" data
79   REAL(wp) ::   rn_pfac     ! multiplication factor for precipitation
80   REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem)
81   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
82   LOGICAL  ::   ln_bulk2z   ! logical flag for case where z(q,t) and z(u) are specified in the namelist
83   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
84   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
85
86   !! * Substitutions
87#  include "domzgr_substitute.h90"
88#  include "vectopt_loop_substitute.h90"
89   !!----------------------------------------------------------------------
90   !! NEMO/OPA 3.3 , NEMO-consortium (2010)
91   !! $Id$
92   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
93   !!----------------------------------------------------------------------
94CONTAINS
95
96   SUBROUTINE sbc_blk_core( kt )
97      !!---------------------------------------------------------------------
98      !!                    ***  ROUTINE sbc_blk_core  ***
99      !!                   
100      !! ** Purpose :   provide at each time step the surface ocean fluxes
101      !!      (momentum, heat, freshwater and runoff)
102      !!
103      !! ** Method  : (1) READ each fluxes in NetCDF files:
104      !!      the 10m wind velocity (i-component) (m/s)    at T-point
105      !!      the 10m wind velocity (j-component) (m/s)    at T-point
106      !!      the 10m or 2m specific humidity     ( % )
107      !!      the solar heat                      (W/m2)
108      !!      the Long wave                       (W/m2)
109      !!      the 10m or 2m air temperature       (Kelvin)
110      !!      the total precipitation (rain+snow) (Kg/m2/s)
111      !!      the snow (solid prcipitation)       (kg/m2/s)
112      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T)
113      !!              (2) CALL blk_oce_core
114      !!
115      !!      C A U T I O N : never mask the surface stress fields
116      !!                      the stress is assumed to be in the (i,j) mesh referential
117      !!
118      !! ** Action  :   defined at each time-step at the air-sea interface
119      !!              - utau, vtau  i- and j-component of the wind stress
120      !!              - taum, wndm  wind stress and 10m wind modules at T-point
121      !!              - qns, qsr    non-solar and solar heat fluxes
122      !!              - emp         upward mass flux (evapo. - precip.)
123      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present)
124      !!                            (set in limsbc(_2).F90)
125      !!----------------------------------------------------------------------
126      INTEGER, INTENT(in) ::   kt   ! ocean time step
127      !!
128      INTEGER  ::   ierror   ! return error code
129      INTEGER  ::   ifpr     ! dummy loop indice
130      INTEGER  ::   jfld     ! dummy loop arguments
131      INTEGER  ::   ios      ! Local integer output status for namelist read
132      !!
133      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files
134      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
135      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
136      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 "
137      TYPE(FLD_N) ::   sn_tdif                                 !   "                                 "
138      NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac, rn_efac, rn_vfac,  &
139         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
140         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
141         &                  sn_tdif, rn_zqt , ln_bulk2z, rn_zu
142      !!---------------------------------------------------------------------
143
144      !                                         ! ====================== !
145      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
146         !                                      ! ====================== !
147         !
148         REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters
149         READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901)
150901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in reference namelist', lwp )
151
152         REWIND( numnam_cfg )              ! Namelist namsbc_core in configuration namelist : CORE bulk parameters
153         READ  ( numnam_cfg, namsbc_core, IOSTAT = ios, ERR = 902 )
154902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in configuration namelist', lwp )
155
156         IF(lwm) WRITE ( numond, namsbc_core )
157         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
158         IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   & 
159            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
160         IF( ln_dm2dc .AND. sn_qsr%ln_tint ) THEN
161            CALL ctl_warn( 'sbc_blk_core: ln_dm2dc is taking care of the temporal interpolation of daily qsr',   &
162                 &         '              ==> We force time interpolation = .false. for qsr' )
163            sn_qsr%ln_tint = .false.
164         ENDIF
165         !                                         ! store namelist information in an array
166         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
167         slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
168         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
169         slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
170         slf_i(jp_tdif) = sn_tdif
171         !                 
172         lhftau = ln_taudif                        ! do we use HF tau information?
173         jfld = jpfld - COUNT( (/.NOT. lhftau/) )
174         !
175         ALLOCATE( sf(jfld), STAT=ierror )         ! set sf structure
176         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core: unable to allocate sf structure' )
177         DO ifpr= 1, jfld
178            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
179            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
180         END DO
181         !                                         ! fill sf with slf_i and control print
182         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
183         !
184         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
185         !
186      ENDIF
187
188      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
189
190      !                                            ! compute the surface ocean fluxes using CORE bulk formulea
191      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( kt, sf, sst_m, ssu_m, ssv_m )
192
193      ! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery
194      IF( ltrcdm2dc )   CALL blk_bio_meanqsr
195
196#if defined key_cice
197      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
198         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
199         qsr_ice(:,:,1)   = sf(jp_qsr)%fnow(:,:,1)
200         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)         
201         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
202         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
203         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
204         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
205         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
206      ENDIF
207#endif
208      !
209   END SUBROUTINE sbc_blk_core
210   
211   
212   SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv )
213      !!---------------------------------------------------------------------
214      !!                     ***  ROUTINE blk_core  ***
215      !!
216      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
217      !!      the ocean surface at each time step
218      !!
219      !! ** Method  :   CORE bulk formulea for the ocean using atmospheric
220      !!      fields read in sbc_read
221      !!
222      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
223      !!              - vtau    : j-component of the stress at V-point  (N/m2)
224      !!              - taum    : Wind stress module at T-point         (N/m2)
225      !!              - wndm    : Wind speed module at T-point          (m/s)
226      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
227      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
228      !!              - evap    : Evaporation over the ocean            (kg/m2/s)
229      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
230      !!
231      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
232      !!---------------------------------------------------------------------
233      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
234      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
235      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
236      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
237      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
238      !
239      INTEGER  ::   ji, jj               ! dummy loop indices
240      REAL(wp) ::   zcoef_qsatw, zztmp   ! local variable
241      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point
242      REAL(wp), DIMENSION(:,:), POINTER ::   zqsatw            ! specific humidity at pst
243      REAL(wp), DIMENSION(:,:), POINTER ::   zqlw, zqsb        ! long wave and sensible heat fluxes
244      REAL(wp), DIMENSION(:,:), POINTER ::   zqla, zevap       ! latent heat fluxes and evaporation
245      REAL(wp), DIMENSION(:,:), POINTER ::   Cd                ! transfer coefficient for momentum      (tau)
246      REAL(wp), DIMENSION(:,:), POINTER ::   Ch                ! transfer coefficient for sensible heat (Q_sens)
247      REAL(wp), DIMENSION(:,:), POINTER ::   Ce                ! tansfert coefficient for evaporation   (Q_lat)
248      REAL(wp), DIMENSION(:,:), POINTER ::   zst               ! surface temperature in Kelvin
249      REAL(wp), DIMENSION(:,:), POINTER ::   zt_zu             ! air temperature at wind speed height
250      REAL(wp), DIMENSION(:,:), POINTER ::   zq_zu             ! air spec. hum.  at wind speed height
251      !!---------------------------------------------------------------------
252      !
253      IF( nn_timing == 1 )  CALL timing_start('blk_oce_core')
254      !
255      CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap )
256      CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
257      !
258      ! local scalars ( place there for vector optimisation purposes)
259      zcoef_qsatw = 0.98 * 640380. / rhoa
260     
261      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
262
263      ! ----------------------------------------------------------------------------- !
264      !      0   Wind components and module at T-point relative to the moving ocean   !
265      ! ----------------------------------------------------------------------------- !
266
267      ! ... components ( U10m - U_oce ) at T-point (unmasked)
268      zwnd_i(:,:) = 0.e0 
269      zwnd_j(:,:) = 0.e0
270#if defined key_cyclone
271# if defined key_vectopt_loop
272!CDIR COLLAPSE
273# endif
274      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add Manu !
275      DO jj = 2, jpjm1
276         DO ji = fs_2, fs_jpim1   ! vect. opt.
277            sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
278            sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
279         END DO
280      END DO
281#endif
282#if defined key_vectopt_loop
283!CDIR COLLAPSE
284#endif
285      DO jj = 2, jpjm1
286         DO ji = fs_2, fs_jpim1   ! vect. opt.
287            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
288            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
289         END DO
290      END DO
291      CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
292      CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
293      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
294!CDIR NOVERRCHK
295!CDIR COLLAPSE
296      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
297         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
298
299      ! ----------------------------------------------------------------------------- !
300      !      I   Radiative FLUXES                                                     !
301      ! ----------------------------------------------------------------------------- !
302   
303      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
304      zztmp = 1. - albo
305      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
306      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
307      ENDIF
308!CDIR COLLAPSE
309      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
310      ! ----------------------------------------------------------------------------- !
311      !     II    Turbulent FLUXES                                                    !
312      ! ----------------------------------------------------------------------------- !
313
314      ! ... specific humidity at SST and IST
315!CDIR NOVERRCHK
316!CDIR COLLAPSE
317      zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 
318
319      ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point :
320      IF( ln_2m ) THEN
321         !! If air temp. and spec. hum. are given at different height (2m) than wind (10m) :
322         CALL TURB_CORE_2Z(2.,10., zst   , sf(jp_tair)%fnow,         &
323            &                      zqsatw, sf(jp_humi)%fnow, wndm,   &
324            &                      Cd    , Ch              , Ce  ,   &
325            &                      zt_zu , zq_zu                   )
326      ELSE IF( ln_bulk2z ) THEN
327         !! If the height of the air temp./spec. hum. and wind are to be specified by hand :
328         IF( rn_zqt == rn_zu ) THEN
329            !! If air temp. and spec. hum. are at the same height as wind :
330            CALL TURB_CORE_1Z( rn_zu, zst   , sf(jp_tair)%fnow(:,:,1),       &
331               &                      zqsatw, sf(jp_humi)%fnow(:,:,1), wndm, &
332               &                      Cd    , Ch                     , Ce  )
333         ELSE
334            !! If air temp. and spec. hum. are at a different height to wind :
335            CALL TURB_CORE_2Z(rn_zqt, rn_zu , zst   , sf(jp_tair)%fnow,         &
336               &                              zqsatw, sf(jp_humi)%fnow, wndm,   &
337               &                              Cd    , Ch              , Ce  ,   &
338               &                              zt_zu , zq_zu                 )
339         ENDIF
340      ELSE
341         !! If air temp. and spec. hum. are given at same height than wind (10m) :
342!gm bug?  at the compiling phase, add a copy in temporary arrays...  ==> check perf
343!         CALL TURB_CORE_1Z( 10., zst   (:,:), sf(jp_tair)%fnow(:,:),              &
344!            &                    zqsatw(:,:), sf(jp_humi)%fnow(:,:), wndm(:,:),   &
345!            &                    Cd    (:,:),             Ch  (:,:), Ce  (:,:)  )
346!gm bug
347! ARPDBG - this won't compile with gfortran. Fix but check performance
348! as per comment above.
349         CALL TURB_CORE_1Z( 10., zst   , sf(jp_tair)%fnow(:,:,1),       &
350            &                    zqsatw, sf(jp_humi)%fnow(:,:,1), wndm, &
351            &                    Cd    , Ch              , Ce    )
352      ENDIF
353
354      ! ... tau module, i and j component
355      DO jj = 1, jpj
356         DO ji = 1, jpi
357            zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj)
358            taum  (ji,jj) = zztmp * wndm  (ji,jj)
359            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
360            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
361         END DO
362      END DO
363
364      ! ... add the HF tau contribution to the wind stress module?
365      IF( lhftau ) THEN 
366!CDIR COLLAPSE
367         taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
368      ENDIF
369      CALL iom_put( "taum_oce", taum )   ! output wind stress module
370
371      ! ... utau, vtau at U- and V_points, resp.
372      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
373      DO jj = 1, jpjm1
374         DO ji = 1, fs_jpim1
375            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) )
376            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) )
377         END DO
378      END DO
379      CALL lbc_lnk( utau(:,:), 'U', -1. )
380      CALL lbc_lnk( vtau(:,:), 'V', -1. )
381
382      !  Turbulent fluxes over ocean
383      ! -----------------------------
384      IF( ln_2m .OR. ( ln_bulk2z .AND. rn_zqt /= rn_zu ) ) THEN
385         ! Values of temp. and hum. adjusted to height of wind must be used
386         zevap(:,:) = rn_efac * MAX( 0.e0, rhoa    *Ce(:,:)*( zqsatw(:,:) - zq_zu(:,:) ) * wndm(:,:) )   ! Evaporation
387         zqsb (:,:) =                      rhoa*cpa*Ch(:,:)*( zst   (:,:) - zt_zu(:,:) ) * wndm(:,:)     ! Sensible Heat
388      ELSE
389!CDIR COLLAPSE
390         zevap(:,:) = rn_efac * MAX( 0.e0, rhoa    *Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) ) * wndm(:,:) )   ! Evaporation
391!CDIR COLLAPSE
392         zqsb (:,:) =            rhoa*cpa*Ch(:,:)*( zst   (:,:) - sf(jp_tair)%fnow(:,:,1) ) * wndm(:,:)     ! Sensible Heat
393      ENDIF
394!CDIR COLLAPSE
395      zqla (:,:) = Lv * zevap(:,:)                                                              ! Latent Heat
396
397      IF(ln_ctl) THEN
398         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce_core: zqla   : ', tab2d_2=Ce , clinfo2=' Ce  : ' )
399         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce_core: zqsb   : ', tab2d_2=Ch , clinfo2=' Ch  : ' )
400         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce_core: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
401         CALL prt_ctl( tab2d_1=zqsatw, clinfo1=' blk_oce_core: zqsatw : ', tab2d_2=zst, clinfo2=' zst : ' )
402         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce_core: utau   : ', mask1=umask,   &
403            &          tab2d_2=vtau  , clinfo2=              ' vtau : '  , mask2=vmask )
404         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_core: wndm   : ')
405         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce_core: zst    : ')
406      ENDIF
407       
408      ! ----------------------------------------------------------------------------- !
409      !     III    Total FLUXES                                                       !
410      ! ----------------------------------------------------------------------------- !
411     
412!CDIR COLLAPSE
413      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
414         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
415!CDIR COLLAPSE
416      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar flux
417         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
418         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
419         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
420         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &   
421         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
422         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic 
423      !
424      CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean
425      CALL iom_put( "qsb_oce", - zqsb )                 ! output downward sensible heat over the ocean
426      CALL iom_put( "qla_oce", - zqla )                 ! output downward latent   heat over the ocean
427      CALL iom_put( "qhc_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
428      CALL iom_put( "qns_oce",   qns  )                 ! output downward non solar heat over the ocean
429      !
430      IF(ln_ctl) THEN
431         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce_core: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
432         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_core: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
433         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce_core: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
434         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce_core: utau   : ', mask1=umask,   &
435            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
436      ENDIF
437      !
438      CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap )
439      CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
440      !
441      IF( nn_timing == 1 )  CALL timing_stop('blk_oce_core')
442      !
443   END SUBROUTINE blk_oce_core
444 
445   SUBROUTINE blk_bio_meanqsr
446      !!---------------------------------------------------------------------
447      !!                     ***  ROUTINE blk_bio_meanqsr
448      !!                     
449      !! ** Purpose :   provide daily qsr_mean for PISCES when
450      !!                analytic diurnal cycle is applied in physic
451      !!               
452      !! ** Method  :   add part where there is no ice
453      !!
454      !!---------------------------------------------------------------------
455      IF( nn_timing == 1 )  CALL timing_start('blk_bio_meanqsr')
456
457      qsr_mean(:,:) = (1. - albo ) *  sf(jp_qsr)%fnow(:,:,1)
458
459      IF( nn_timing == 1 )  CALL timing_stop('blk_bio_meanqsr')
460
461   END SUBROUTINE blk_bio_meanqsr
462 
463 
464   SUBROUTINE blk_ice_meanqsr(palb,p_qsr_mean,pdim)
465      !!---------------------------------------------------------------------
466      !!
467      !! ** Purpose :   provide the daily qsr_mean over sea_ice for PISCES when
468      !!                analytic diurnal cycle is applied in physic
469      !!
470      !! ** Method  :   compute qsr
471      !!
472      !!---------------------------------------------------------------------
473      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb       ! ice albedo (clear sky) (alb_ice_cs)               [%]
474      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr_mean !     solar heat flux over ice (T-point)         [W/m2]
475      INTEGER                   , INTENT(in   ) ::   pdim       ! number of ice categories
476      !!
477      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays)
478      INTEGER  ::   ji, jj, jl    ! dummy loop indices
479      REAL(wp) ::   zztmp         ! temporary variable
480      !!---------------------------------------------------------------------
481      IF( nn_timing == 1 )  CALL timing_start('blk_ice_meanqsr')
482      !
483      ijpl  = pdim                            ! number of ice categories
484      zztmp = 1. / ( 1. - albo )
485      !                                     ! ========================== !
486      DO jl = 1, ijpl                       !  Loop over ice categories  !
487         !                                  ! ========================== !
488         DO jj = 1 , jpj
489            DO ji = 1, jpi
490                  p_qsr_mean(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr_mean(ji,jj)
491            END DO
492         END DO
493      END DO
494      !
495      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_meanqsr')
496      !
497   END SUBROUTINE blk_ice_meanqsr 
498 
499   
500   SUBROUTINE blk_ice_core(  pst   , pui   , pvi   , palb ,   &
501      &                      p_taui, p_tauj, p_qns , p_qsr,   &
502      &                      p_qla , p_dqns, p_dqla,          &
503      &                      p_tpr , p_spr ,                  &
504      &                      p_fr1 , p_fr2 , cd_grid, pdim  ) 
505      !!---------------------------------------------------------------------
506      !!                     ***  ROUTINE blk_ice_core  ***
507      !!
508      !! ** Purpose :   provide the surface boundary condition over sea-ice
509      !!
510      !! ** Method  :   compute momentum, heat and freshwater exchanged
511      !!                between atmosphere and sea-ice using CORE bulk
512      !!                formulea, ice variables and read atmmospheric fields.
513      !!                NB: ice drag coefficient is assumed to be a constant
514      !!
515      !! caution : the net upward water flux has with mm/day unit
516      !!---------------------------------------------------------------------
517      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin]
518      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pui      ! ice surface velocity (i- and i- components      [m/s]
519      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid)
520      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%]
521      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2]
522      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid)
523      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2]
524      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2]
525      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2]
526      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2]
527      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2]
528      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s]
529      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s]
530      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%]
531      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%]
532      CHARACTER(len=1)          , INTENT(in   ) ::   cd_grid  ! ice grid ( C or B-grid)
533      INTEGER                   , INTENT(in   ) ::   pdim     ! number of ice categories
534      !!
535      INTEGER  ::   ji, jj, jl    ! dummy loop indices
536      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays)
537      REAL(wp) ::   zst2, zst3
538      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb
539      REAL(wp) ::   zztmp                                        ! temporary variable
540      REAL(wp) ::   zcoef_frca                                   ! fractional cloud amount
541      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                  ! relative wind module and components at F-point
542      REAL(wp) ::             zwndi_t , zwndj_t                  ! relative wind components at T-point
543      !!
544      REAL(wp), DIMENSION(:,:)  , POINTER ::   z_wnds_t          ! wind speed ( = | U10m - U_ice | ) at T-point
545      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice
546      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice
547      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice
548      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice
549      !!---------------------------------------------------------------------
550      !
551      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core')
552      !
553      CALL wrk_alloc( jpi,jpj, z_wnds_t )
554      CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
555
556      ijpl  = pdim                            ! number of ice categories
557
558      ! local scalars ( place there for vector optimisation purposes)
559      zcoef_wnorm  = rhoa * Cice
560      zcoef_wnorm2 = rhoa * Cice * 0.5
561      zcoef_dqlw   = 4.0 * 0.95 * Stef
562      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8)
563      zcoef_dqsb   = rhoa * cpa * Cice
564      zcoef_frca   = 1.0  - 0.3
565
566!!gm brutal....
567      z_wnds_t(:,:) = 0.e0
568      p_taui  (:,:) = 0.e0
569      p_tauj  (:,:) = 0.e0
570!!gm end
571
572#if defined key_lim3
573      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)   ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init
574#endif
575      ! ----------------------------------------------------------------------------- !
576      !    Wind components and module relative to the moving ocean ( U10m - U_ice )   !
577      ! ----------------------------------------------------------------------------- !
578      SELECT CASE( cd_grid )
579      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
580         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked)
581!CDIR NOVERRCHK
582         DO jj = 2, jpjm1
583            DO ji = 2, jpim1   ! B grid : NO vector opt
584               ! ... scalar wind at I-point (fld being at T-point)
585               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   &
586                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * pui(ji,jj)
587               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   &
588                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * pvi(ji,jj)
589               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
590               ! ... ice stress at I-point
591               p_taui(ji,jj) = zwnorm_f * zwndi_f
592               p_tauj(ji,jj) = zwnorm_f * zwndj_f
593               ! ... scalar wind at T-point (fld being at T-point)
594               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pui(ji,jj+1) + pui(ji+1,jj+1)   &
595                  &                                                    + pui(ji,jj  ) + pui(ji+1,jj  )  )
596               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pvi(ji,jj+1) + pvi(ji+1,jj+1)   &
597                  &                                                    + pvi(ji,jj  ) + pvi(ji+1,jj  )  )
598               z_wnds_t(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
599            END DO
600         END DO
601         CALL lbc_lnk( p_taui  , 'I', -1. )
602         CALL lbc_lnk( p_tauj  , 'I', -1. )
603         CALL lbc_lnk( z_wnds_t, 'T',  1. )
604         !
605      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean)
606#if defined key_vectopt_loop
607!CDIR COLLAPSE
608#endif
609         DO jj = 2, jpj
610            DO ji = fs_2, jpi   ! vect. opt.
611               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pui(ji-1,jj  ) + pui(ji,jj) )  )
612               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pvi(ji  ,jj-1) + pvi(ji,jj) )  )
613               z_wnds_t(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
614            END DO
615         END DO
616#if defined key_vectopt_loop
617!CDIR COLLAPSE
618#endif
619         DO jj = 2, jpjm1
620            DO ji = fs_2, fs_jpim1   ! vect. opt.
621               p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj  ) + z_wnds_t(ji,jj) )                          &
622                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) )
623               p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1  ) + z_wnds_t(ji,jj) )                          &
624                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * pvi(ji,jj) )
625            END DO
626         END DO
627         CALL lbc_lnk( p_taui  , 'U', -1. )
628         CALL lbc_lnk( p_tauj  , 'V', -1. )
629         CALL lbc_lnk( z_wnds_t, 'T',  1. )
630         !
631      END SELECT
632
633      zztmp = 1. / ( 1. - albo )
634      !                                     ! ========================== !
635      DO jl = 1, ijpl                       !  Loop over ice categories  !
636         !                                  ! ========================== !
637!CDIR NOVERRCHK
638!CDIR COLLAPSE
639         DO jj = 1 , jpj
640!CDIR NOVERRCHK
641            DO ji = 1, jpi
642               ! ----------------------------!
643               !      I   Radiative FLUXES   !
644               ! ----------------------------!
645               zst2 = pst(ji,jj,jl) * pst(ji,jj,jl)
646               zst3 = pst(ji,jj,jl) * zst2
647               ! Short Wave (sw)
648               p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
649               ! Long  Wave (lw)
650               ! iovino
651               IF( ff(ji,jj) .GT. 0._wp ) THEN
652                  z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
653               ELSE
654                  z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
655               ENDIF
656               ! lw sensitivity
657               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                               
658
659               ! ----------------------------!
660               !     II    Turbulent FLUXES  !
661               ! ----------------------------!
662
663               ! ... turbulent heat fluxes
664               ! Sensible Heat
665               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * z_wnds_t(ji,jj) * ( pst(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
666               ! Latent Heat
667               p_qla(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * z_wnds_t(ji,jj)   &                           
668                  &                         * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
669               ! Latent heat sensitivity for ice (Dqla/Dt)
670               p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) )
671               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
672               z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj)
673
674               ! ----------------------------!
675               !     III    Total FLUXES     !
676               ! ----------------------------!
677               ! Downward Non Solar flux
678               p_qns (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla (ji,jj,jl)     
679               ! Total non solar heat flux sensitivity for ice
680               p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) )   
681            END DO
682            !
683         END DO
684         !
685      END DO
686      !
687      !--------------------------------------------------------------------
688      ! FRACTIONs of net shortwave radiation which is not absorbed in the
689      ! thin surface layer and penetrates inside the ice cover
690      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
691   
692!CDIR COLLAPSE
693      p_fr1(:,:) = ( 0.18 * ( 1.0 - zcoef_frca ) + 0.35 * zcoef_frca )
694!CDIR COLLAPSE
695      p_fr2(:,:) = ( 0.82 * ( 1.0 - zcoef_frca ) + 0.65 * zcoef_frca )
696       
697!CDIR COLLAPSE
698      p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s]
699!CDIR COLLAPSE
700      p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s]
701      CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation
702      CALL iom_put( 'precip', p_tpr * 86400. )                   ! Total precipitation
703      !
704      IF(ln_ctl) THEN
705         CALL prt_ctl(tab3d_1=p_qla   , clinfo1=' blk_ice_core: p_qla  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb  : ', kdim=ijpl)
706         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw  : ', tab3d_2=p_dqla  , clinfo2=' p_dqla : ', kdim=ijpl)
707         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw : ', kdim=ijpl)
708         CALL prt_ctl(tab3d_1=p_dqns  , clinfo1=' blk_ice_core: p_dqns : ', tab3d_2=p_qsr   , clinfo2=' p_qsr  : ', kdim=ijpl)
709         CALL prt_ctl(tab3d_1=pst     , clinfo1=' blk_ice_core: pst    : ', tab3d_2=p_qns   , clinfo2=' p_qns  : ', kdim=ijpl)
710         CALL prt_ctl(tab2d_1=p_tpr   , clinfo1=' blk_ice_core: p_tpr  : ', tab2d_2=p_spr   , clinfo2=' p_spr  : ')
711         CALL prt_ctl(tab2d_1=p_taui  , clinfo1=' blk_ice_core: p_taui : ', tab2d_2=p_tauj  , clinfo2=' p_tauj : ')
712         CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ')
713      ENDIF
714
715      CALL wrk_dealloc( jpi,jpj, z_wnds_t )
716      CALL wrk_dealloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
717      !
718      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core')
719      !
720   END SUBROUTINE blk_ice_core
721 
722
723   SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a,   &
724      &                        dU , Cd , Ch   , Ce   )
725      !!----------------------------------------------------------------------
726      !!                      ***  ROUTINE  turb_core  ***
727      !!
728      !! ** Purpose :   Computes turbulent transfert coefficients of surface
729      !!                fluxes according to Large & Yeager (2004)
730      !!
731      !! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D
732      !!      Momentum, Latent and sensible heat exchange coefficients
733      !!      Caution: this procedure should only be used in cases when air
734      !!      temperature (T_air), air specific humidity (q_air) and wind (dU)
735      !!      are provided at the same height 'zzu'!
736      !!
737      !! References :   Large & Yeager, 2004 : ???
738      !!----------------------------------------------------------------------
739      REAL(wp)                , INTENT(in   ) ::   zu      ! altitude of wind measurement       [m]
740      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   sst     ! sea surface temperature         [Kelvin]
741      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   T_a     ! potential air temperature       [Kelvin]
742      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   q_sat   ! sea surface specific humidity   [kg/kg]
743      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   q_a     ! specific air humidity           [kg/kg]
744      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   dU      ! wind module |U(zu)-U(0)|        [m/s]
745      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   Cd      ! transfert coefficient for momentum       (tau)
746      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   Ch      ! transfert coefficient for temperature (Q_sens)
747      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   Ce      ! transfert coefficient for evaporation  (Q_lat)
748      !!
749      INTEGER :: j_itt
750      INTEGER , PARAMETER ::   nb_itt = 3
751      REAL(wp), PARAMETER ::   grav   = 9.8   ! gravity                       
752      REAL(wp), PARAMETER ::   kappa  = 0.4   ! von Karman s constant
753
754      REAL(wp), DIMENSION(:,:), POINTER  ::   dU10          ! dU                                   [m/s]
755      REAL(wp), DIMENSION(:,:), POINTER  ::   dT            ! air/sea temperature differeence      [K]
756      REAL(wp), DIMENSION(:,:), POINTER  ::   dq            ! air/sea humidity difference          [K]
757      REAL(wp), DIMENSION(:,:), POINTER  ::   Cd_n10        ! 10m neutral drag coefficient
758      REAL(wp), DIMENSION(:,:), POINTER  ::   Ce_n10        ! 10m neutral latent coefficient
759      REAL(wp), DIMENSION(:,:), POINTER  ::   Ch_n10        ! 10m neutral sensible coefficient
760      REAL(wp), DIMENSION(:,:), POINTER  ::   sqrt_Cd_n10   ! root square of Cd_n10
761      REAL(wp), DIMENSION(:,:), POINTER  ::   sqrt_Cd       ! root square of Cd
762      REAL(wp), DIMENSION(:,:), POINTER  ::   T_vpot        ! virtual potential temperature        [K]
763      REAL(wp), DIMENSION(:,:), POINTER  ::   T_star        ! turbulent scale of tem. fluct.
764      REAL(wp), DIMENSION(:,:), POINTER  ::   q_star        ! turbulent humidity of temp. fluct.
765      REAL(wp), DIMENSION(:,:), POINTER  ::   U_star        ! turb. scale of velocity fluct.
766      REAL(wp), DIMENSION(:,:), POINTER  ::   L             ! Monin-Obukov length                  [m]
767      REAL(wp), DIMENSION(:,:), POINTER  ::   zeta          ! stability parameter at height zu
768      REAL(wp), DIMENSION(:,:), POINTER  ::   U_n10         ! neutral wind velocity at 10m         [m]   
769      REAL(wp), DIMENSION(:,:), POINTER  ::   xlogt, xct, zpsi_h, zpsi_m
770     
771      INTEGER , DIMENSION(:,:), POINTER  ::   stab          ! 1st guess stability test integer
772      !!----------------------------------------------------------------------
773      !
774      IF( nn_timing == 1 )  CALL timing_start('TURB_CORE_1Z')
775      !
776      CALL wrk_alloc( jpi,jpj, stab )   ! integer
777      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
778      CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m )
779
780      !! * Start
781      !! Air/sea differences
782      dU10 = max(0.5, dU)     ! we don't want to fall under 0.5 m/s
783      dT = T_a - sst          ! assuming that T_a is allready the potential temp. at zzu
784      dq = q_a - q_sat
785      !!   
786      !! Virtual potential temperature
787      T_vpot = T_a*(1. + 0.608*q_a)
788      !!
789      !! Neutral Drag Coefficient
790      stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0
791      IF  ( ln_cdgw ) THEN
792        cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
793        Cd_n10(:,:) =   cdn_wave
794      ELSE
795        Cd_n10  = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a)
796      ENDIF
797      sqrt_Cd_n10 = sqrt(Cd_n10)
798      Ce_n10  = 1.e-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b)
799      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) !   L & Y eq. (6c), (6d)
800      !!
801      !! Initializing transfert coefficients with their first guess neutral equivalents :
802      Cd = Cd_n10 ;  Ce = Ce_n10 ;  Ch = Ch_n10 ;  sqrt_Cd = sqrt(Cd)
803
804      !! * Now starting iteration loop
805      DO j_itt=1, nb_itt
806         !! Turbulent scales :
807         U_star  = sqrt_Cd*dU10                !   L & Y eq. (7a)
808         T_star  = Ch/sqrt_Cd*dT               !   L & Y eq. (7b)
809         q_star  = Ce/sqrt_Cd*dq               !   L & Y eq. (7c)
810
811         !! Estimate the Monin-Obukov length :
812         L  = (U_star**2)/( kappa*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) )
813
814         !! Stability parameters :
815         zeta  = zu/L ;  zeta   = sign( min(abs(zeta),10.0), zeta )
816         zpsi_h  = psi_h(zeta)
817         zpsi_m  = psi_m(zeta)
818
819         IF  ( ln_cdgw ) THEN
820           sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd;
821         ELSE
822           !! Shifting the wind speed to 10m and neutral stability :
823           U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) !  L & Y eq. (9a)
824
825           !! Updating the neutral 10m transfer coefficients :
826           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a)
827           sqrt_Cd_n10 = sqrt(Cd_n10)
828           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b)
829           stab    = 0.5 + sign(0.5,zeta)
830           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab))           !  L & Y eq. (6c), (6d)
831
832           !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) :
833           !!
834           xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m)
835           Cd  = Cd_n10/(xct*xct) ;  sqrt_Cd = sqrt(Cd)
836         ENDIF
837         !!
838         xlogt = log(zu/10.) - zpsi_h
839         !!
840         xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10
841         Ch  = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct
842         !!
843         xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10
844         Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
845         !!
846      END DO
847      !!
848      CALL wrk_dealloc( jpi,jpj, stab )   ! integer
849      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
850      CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m )
851      !
852      IF( nn_timing == 1 )  CALL timing_stop('TURB_CORE_1Z')
853      !
854    END SUBROUTINE TURB_CORE_1Z
855
856
857    SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu)
858      !!----------------------------------------------------------------------
859      !!                      ***  ROUTINE  turb_core  ***
860      !!
861      !! ** Purpose :   Computes turbulent transfert coefficients of surface
862      !!                fluxes according to Large & Yeager (2004).
863      !!
864      !! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D
865      !!      Momentum, Latent and sensible heat exchange coefficients
866      !!      Caution: this procedure should only be used in cases when air
867      !!      temperature (T_air) and air specific humidity (q_air) are at a
868      !!      different height to wind (dU).
869      !!
870      !! References :   Large & Yeager, 2004 : ???
871      !!----------------------------------------------------------------------
872      REAL(wp), INTENT(in   )                     ::   zt       ! height for T_zt and q_zt                   [m]
873      REAL(wp), INTENT(in   )                     ::   zu       ! height for dU                              [m]
874      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature              [Kelvin]
875      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   T_zt     ! potential air temperature            [Kelvin]
876      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_sat    ! sea surface specific humidity         [kg/kg]
877      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                 [kg/kg]
878      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   dU       ! relative wind module |U(zu)-U(0)|       [m/s]
879      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
880      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
881      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
882      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   T_zu     ! air temp. shifted at zu                     [K]
883      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. hum.  shifted at zu               [kg/kg]
884
885      INTEGER :: j_itt
886      INTEGER , PARAMETER :: nb_itt = 5              ! number of itterations
887      REAL(wp), PARAMETER ::   grav   = 9.8          ! gravity                       
888      REAL(wp), PARAMETER ::   kappa  = 0.4          ! von Karman's constant
889     
890      REAL(wp), DIMENSION(:,:), POINTER ::   dU10          ! dU                                [m/s]
891      REAL(wp), DIMENSION(:,:), POINTER ::   dT            ! air/sea temperature differeence   [K]
892      REAL(wp), DIMENSION(:,:), POINTER ::   dq            ! air/sea humidity difference       [K]
893      REAL(wp), DIMENSION(:,:), POINTER ::   Cd_n10        ! 10m neutral drag coefficient
894      REAL(wp), DIMENSION(:,:), POINTER ::   Ce_n10        ! 10m neutral latent coefficient
895      REAL(wp), DIMENSION(:,:), POINTER ::   Ch_n10        ! 10m neutral sensible coefficient
896      REAL(wp), DIMENSION(:,:), POINTER ::   sqrt_Cd_n10   ! root square of Cd_n10
897      REAL(wp), DIMENSION(:,:), POINTER ::   sqrt_Cd       ! root square of Cd
898      REAL(wp), DIMENSION(:,:), POINTER ::   T_vpot        ! virtual potential temperature        [K]
899      REAL(wp), DIMENSION(:,:), POINTER ::   T_star        ! turbulent scale of tem. fluct.
900      REAL(wp), DIMENSION(:,:), POINTER ::   q_star        ! turbulent humidity of temp. fluct.
901      REAL(wp), DIMENSION(:,:), POINTER ::   U_star        ! turb. scale of velocity fluct.
902      REAL(wp), DIMENSION(:,:), POINTER ::   L             ! Monin-Obukov length                  [m]
903      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_u        ! stability parameter at height zu
904      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_t        ! stability parameter at height zt
905      REAL(wp), DIMENSION(:,:), POINTER ::   U_n10         ! neutral wind velocity at 10m        [m]
906      REAL(wp), DIMENSION(:,:), POINTER ::   xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m
907
908      INTEGER , DIMENSION(:,:), POINTER ::   stab          ! 1st stability test integer
909      !!----------------------------------------------------------------------
910      !
911      IF( nn_timing == 1 )  CALL timing_start('TURB_CORE_2Z')
912      !
913      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
914      CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 )
915      CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
916      CALL wrk_alloc( jpi,jpj, stab )   ! interger
917
918      !! Initial air/sea differences
919      dU10 = max(0.5, dU)      !  we don't want to fall under 0.5 m/s
920      dT = T_zt - sst 
921      dq = q_zt - q_sat
922
923      !! Neutral Drag Coefficient :
924      stab = 0.5 + sign(0.5,dT)                 ! stab = 1  if dT > 0  -> STABLE
925      IF( ln_cdgw ) THEN
926        cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
927        Cd_n10(:,:) =   cdn_wave
928      ELSE
929        Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 
930      ENDIF
931      sqrt_Cd_n10 = sqrt(Cd_n10)
932      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 )
933      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))
934
935      !! Initializing transf. coeff. with their first guess neutral equivalents :
936      Cd = Cd_n10 ;  Ce = Ce_n10 ;  Ch = Ch_n10 ;  sqrt_Cd = sqrt(Cd)
937
938      !! Initializing z_u values with z_t values :
939      T_zu = T_zt ;  q_zu = q_zt
940
941      !!  * Now starting iteration loop
942      DO j_itt=1, nb_itt
943         dT = T_zu - sst ;  dq = q_zu - q_sat ! Updating air/sea differences
944         T_vpot = T_zu*(1. + 0.608*q_zu)    ! Updating virtual potential temperature at zu
945         U_star = sqrt_Cd*dU10                ! Updating turbulent scales :   (L & Y eq. (7))
946         T_star  = Ch/sqrt_Cd*dT              !
947         q_star  = Ce/sqrt_Cd*dq              !
948         !!
949         L = (U_star*U_star) &                ! Estimate the Monin-Obukov length at height zu
950              & / (kappa*grav/T_vpot*(T_star*(1.+0.608*q_zu) + 0.608*T_zu*q_star))
951         !! Stability parameters :
952         zeta_u  = zu/L  ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u )
953         zeta_t  = zt/L  ;  zeta_t = sign( min(abs(zeta_t),10.0), zeta_t )
954         zpsi_hu = psi_h(zeta_u)
955         zpsi_ht = psi_h(zeta_t)
956         zpsi_m  = psi_m(zeta_u)
957         !!
958         !! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a))
959!        U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)))
960         U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m))
961         !!
962         !! Shifting temperature and humidity at zu :          (L & Y eq. (9b-9c))
963!        T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t))
964         T_zu = T_zt - T_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht)
965!        q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t))
966         q_zu = q_zt - q_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht)
967         !!
968         !! q_zu cannot have a negative value : forcing 0
969         stab = 0.5 + sign(0.5,q_zu) ;  q_zu = stab*q_zu
970         !!
971         IF( ln_cdgw ) THEN
972            sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd;
973         ELSE
974           !! Updating the neutral 10m transfer coefficients :
975           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a)
976           sqrt_Cd_n10 = sqrt(Cd_n10)
977           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b)
978           stab    = 0.5 + sign(0.5,zeta_u)
979           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) ! L & Y eq. (6c-6d)
980           !!
981           !!
982           !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) :
983           xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)   ! L & Y eq. (10a)
984           Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd)
985         ENDIF
986         !!
987         xlogt = log(zu/10.) - zpsi_hu
988         !!
989         xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10               ! L & Y eq. (10b)
990         Ch  = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct
991         !!
992         xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10               ! L & Y eq. (10c)
993         Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
994         !!
995         !!
996      END DO
997      !!
998      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
999      CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 )
1000      CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
1001      CALL wrk_dealloc( jpi,jpj, stab )   ! interger
1002      !
1003      IF( nn_timing == 1 )  CALL timing_stop('TURB_CORE_2Z')
1004      !
1005    END SUBROUTINE TURB_CORE_2Z
1006
1007
1008    FUNCTION psi_m(zta)   !! Psis, L & Y eq. (8c), (8d), (8e)
1009      !-------------------------------------------------------------------------------
1010      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta
1011
1012      REAL(wp), PARAMETER :: pi = 3.141592653589793_wp
1013      REAL(wp), DIMENSION(jpi,jpj)             :: psi_m
1014      REAL(wp), DIMENSION(:,:), POINTER        :: X2, X, stabit
1015      !-------------------------------------------------------------------------------
1016
1017      CALL wrk_alloc( jpi,jpj, X2, X, stabit )
1018
1019      X2 = sqrt(abs(1. - 16.*zta))  ;  X2 = max(X2 , 1.0) ;  X  = sqrt(X2)
1020      stabit    = 0.5 + sign(0.5,zta)
1021      psi_m = -5.*zta*stabit  &                                                          ! Stable
1022         &    + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2)  ! Unstable
1023
1024      CALL wrk_dealloc( jpi,jpj, X2, X, stabit )
1025      !
1026    END FUNCTION psi_m
1027
1028
1029    FUNCTION psi_h( zta )    !! Psis, L & Y eq. (8c), (8d), (8e)
1030      !-------------------------------------------------------------------------------
1031      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zta
1032      !
1033      REAL(wp), DIMENSION(jpi,jpj)             ::   psi_h
1034      REAL(wp), DIMENSION(:,:), POINTER        :: X2, X, stabit
1035      !-------------------------------------------------------------------------------
1036
1037      CALL wrk_alloc( jpi,jpj, X2, X, stabit )
1038
1039      X2 = sqrt(abs(1. - 16.*zta))  ;  X2 = max(X2 , 1.) ;  X  = sqrt(X2)
1040      stabit    = 0.5 + sign(0.5,zta)
1041      psi_h = -5.*zta*stabit  &                                       ! Stable
1042         &    + (1. - stabit)*(2.*log( (1. + X2)/2. ))                 ! Unstable
1043
1044      CALL wrk_dealloc( jpi,jpj, X2, X, stabit )
1045      !
1046    END FUNCTION psi_h
1047 
1048   !!======================================================================
1049END MODULE sbcblk_core
Note: See TracBrowser for help on using the repository browser.