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

source: branches/2011/dev_r2802_UKMO8_cice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 @ 2874

Last change on this file since 2874 was 2874, checked in by charris, 13 years ago

Code for running NEMO with CICE (for fully coupled mode this should be used in combination with dev_r2802_UKMO8_sbccpl). Changes are described briefly below.

physct: Constants modified to be consistent with CICE

nemogcm / prtctl / mppini: Changes to NEMO decomposition (activated using key_nemocice_decomp) to produce 'square' options in CICE. Can run without this key / code but this requires a global gather / scatter in the NEMO-CICE coupling which gets very slow on large processors numbers.

sbc_ice: CICE options and arrays added

sbcmod: CICE option added, including calls for initialising and finalising CICE.

sbcblk_core: Make sure necessary forcing field are available for CICE

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