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 trunk/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMO/OPA_SRC/SBC/sbcblk_core.F90 @ 1200

Last change on this file since 1200 was 1200, checked in by rblod, 16 years ago

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