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

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 @ 8879

Last change on this file since 8879 was 8879, checked in by frrh, 6 years ago

Merge in http://fcm3/projects/NEMO.xm/log/branches/UKMO/dev_r8183_ICEMODEL_svn_removed
revisions 8738:8847 inclusive.

File size: 73.0 KB
Line 
1MODULE sbcblk
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!                         Aerodynamic Bulk Formulas
6   !!                        SUCCESSOR OF "sbcblk_core"
7   !!=====================================================================
8   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original CORE code
9   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier)  improved CORE bulk and its user interface
10   !!            3.0  !  2006-06  (G. Madec)  sbc rewritting
11   !!             -   !  2006-12  (L. Brodeau)  Original code for turb_core
12   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
13   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
14   !!            3.4  !  2011-11  (C. Harris)  Fill arrays required by CICE
15   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
16   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
17   !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
18   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle)
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
24   !!   sbc_blk       : bulk formulation as ocean surface boundary condition
25   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
26   !!   blk_ice       : computes momentum, heat and freshwater fluxes over sea ice
27   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
28   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
29   !!   q_sat         : saturation humidity as a function of SLP and temperature
30   !!   L_vap         : latent heat of vaporization of water as a function of temperature
31   !!----------------------------------------------------------------------
32   USE oce            ! ocean dynamics and tracers
33   USE dom_oce        ! ocean space and time domain
34   USE phycst         ! physical constants
35   USE fldread        ! read input fields
36   USE sbc_oce        ! Surface boundary condition: ocean fields
37   USE cyclone        ! Cyclone 10m wind form trac of cyclone centres
38   USE sbcdcy         ! surface boundary condition: diurnal cycle
39   USE sbcwave , ONLY :   cdn_wave ! wave module
40   USE sbc_ice        ! Surface boundary condition: ice fields
41   USE lib_fortran    ! to use key_nosignedzero
42#if defined key_lim3
43   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su
44   USE icethd_dh      ! for CALL ice_thd_snwblow
45   USE icethd_zdf, ONLY:   rn_cnd_s  ! for blk_ice_qcn
46#endif
47   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
48   USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003)
49   USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)
50   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31)
51   !
52   USE iom            ! I/O manager library
53   USE in_out_manager ! I/O manager
54   USE lib_mpp        ! distribued memory computing library
55   USE wrk_nemo       ! work arrays
56   USE timing         ! Timing
57   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
58   USE prtctl         ! Print control
59
60   IMPLICIT NONE
61   PRIVATE
62
63   PUBLIC   sbc_blk_init  ! called in sbcmod
64   PUBLIC   sbc_blk       ! called in sbcmod
65#if defined key_lim3
66   PUBLIC   blk_ice_tau   ! routine called in icestp module
67   PUBLIC   blk_ice_flx   ! routine called in icestp module
68   PUBLIC   blk_ice_qcn   ! routine called in icestp module
69#endif 
70
71!!Lolo: should ultimately be moved in the module with all physical constants ?
72!!gm  : In principle, yes.
73   REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg]
74   REAL(wp), PARAMETER ::   Cp_vap = 1860.0       !: Specic heat of water vapor, constant pressure  [J/K/kg]
75   REAL(wp), PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg]
76   REAL(wp), PARAMETER ::   R_vap = 461.495_wp    !: Specific gas constant for water vapor          [J/K/kg]
77   REAL(wp), PARAMETER ::   reps0 = R_dry/R_vap   !: ratio of gas constant for dry air and water vapor => ~ 0.622
78   REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
79
80   INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read
81   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
82   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
83   INTEGER , PARAMETER ::   jp_tair = 3           ! index of 10m air temperature             (Kelvin)
84   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % )
85   INTEGER , PARAMETER ::   jp_qsr  = 5           ! index of solar heat                      (W/m2)
86   INTEGER , PARAMETER ::   jp_qlw  = 6           ! index of Long wave                       (W/m2)
87   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
88   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
89   INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa)
90   INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point
91
92   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
93
94   !                                             !!! Bulk parameters
95   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air (only used for ice fluxes now...)
96   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation
97   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant
98   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice
99   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant
100   !
101   !                           !!* Namelist namsbc_blk : bulk parameters
102   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
103   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
104   LOGICAL  ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013)
105   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
106   !
107   LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data
108   REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation
109   REAL(wp) ::   rn_efac        ! multiplication factor for evaporation (clem)
110   REAL(wp) ::   rn_vfac        ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
111   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements
112   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements
113   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
114   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
115   !
116   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau)
117   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens)
118   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat)
119   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
120   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme)
121   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
122
123   INTEGER  ::   nblk           ! choice of the bulk algorithm
124   !                            ! associated indices:
125   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008)
126   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
127   INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013)
128   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31)
129
130   !! * Substitutions
131#  include "vectopt_loop_substitute.h90"
132   !!----------------------------------------------------------------------
133   !! NEMO/OPA 3.7 , NEMO-consortium (2014)
134   !! $Id$
135   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
136   !!----------------------------------------------------------------------
137CONTAINS
138
139   INTEGER FUNCTION sbc_blk_alloc()
140      !!-------------------------------------------------------------------
141      !!             ***  ROUTINE sbc_blk_alloc ***
142      !!-------------------------------------------------------------------
143      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
144         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
145      !
146      IF( lk_mpp             )   CALL mpp_sum ( sbc_blk_alloc )
147      IF( sbc_blk_alloc /= 0 )   CALL ctl_warn('sbc_blk_alloc: failed to allocate arrays')
148   END FUNCTION sbc_blk_alloc
149
150   SUBROUTINE sbc_blk_init
151      !!---------------------------------------------------------------------
152      !!                    ***  ROUTINE sbc_blk_init  ***
153      !!
154      !! ** Purpose :   choose and initialize a bulk formulae formulation
155      !!
156      !! ** Method  :
157      !!
158      !!      C A U T I O N : never mask the surface stress fields
159      !!                      the stress is assumed to be in the (i,j) mesh referential
160      !!
161      !! ** Action  :   
162      !!
163      !!----------------------------------------------------------------------
164      INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument
165      INTEGER  ::   ios, ierror, ioptio   ! Local integer
166      !!
167      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
168      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
169      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
170      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
171      TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        "
172      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
173         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
174         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
175         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
176         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
177      !!---------------------------------------------------------------------
178      !
179      !                                      ! allocate sbc_blk_core array
180      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
181      !
182      !                             !** read bulk namelist 
183      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters
184      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
185901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp )
186      !
187      REWIND( numnam_cfg )                !* Namelist namsbc_blk in configuration namelist : bulk parameters
188      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
189902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist', lwp )
190      !
191      IF(lwm) WRITE( numond, namsbc_blk )
192      !
193      !                             !** initialization of the chosen bulk formulae (+ check)
194      !                                   !* select the bulk chosen in the namelist and check the choice
195      ;                                                        ioptio = 0
196      IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF
197      IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF
198      IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF
199      IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF
200      !
201      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
202      !
203      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
204         IF( sn_qsr%nfreqh /= 24 )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
205         IF( sn_qsr%ln_tint ) THEN
206            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   &
207               &           '              ==> We force time interpolation = .false. for qsr' )
208            sn_qsr%ln_tint = .false.
209         ENDIF
210      ENDIF
211      !                                   !* set the bulk structure
212      !                                      !- store namelist information in an array
213      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
214      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
215      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
216      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
217      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif
218      !
219      lhftau = ln_taudif                     !- add an extra field if HF stress is used
220      jfld = jpfld - COUNT( (/.NOT.lhftau/) )
221      !
222      !                                      !- allocate the bulk structure
223      ALLOCATE( sf(jfld), STAT=ierror )
224      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
225      DO ifpr= 1, jfld
226         ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
227         IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
228      END DO
229      !                                      !- fill the bulk structure with namelist informations
230      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
231      !
232      IF ( ln_wave ) THEN
233      !Activated wave module but neither drag nor stokes drift activated
234         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) )   THEN
235            CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F')
236      !drag coefficient read from wave model definable only with mfs bulk formulae and core
237         ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN       
238             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
239         ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN
240             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
241         ENDIF
242      ELSE
243      IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor )                & 
244         &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
245         &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
246         &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
247         &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      & 
248         &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
249      ENDIF 
250      !
251      !           
252      IF(lwp) THEN                     !** Control print
253         !
254         WRITE(numout,*)                  !* namelist
255         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):'
256         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
257         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
258         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0
259         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
260         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
261         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt
262         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
263         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
264         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
265         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
266         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
267         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12
268         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15
269         !
270         WRITE(numout,*)
271         SELECT CASE( nblk )              !* Print the choice of bulk algorithm
272         CASE( np_NCAR      )   ;   WRITE(numout,*) '      ===>>   "NCAR" algorithm        (Large and Yeager 2008)'
273         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '      ===>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
274         CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '      ===>>   "COARE 3.5" algorithm   (Edson et al. 2013)'
275         CASE( np_ECMWF     )   ;   WRITE(numout,*) '      ===>>   "ECMWF" algorithm       (IFS cycle 31)'
276         END SELECT
277         !
278      ENDIF
279      !
280   END SUBROUTINE sbc_blk_init
281
282
283   SUBROUTINE sbc_blk( kt )
284      !!---------------------------------------------------------------------
285      !!                    ***  ROUTINE sbc_blk  ***
286      !!
287      !! ** Purpose :   provide at each time step the surface ocean fluxes
288      !!      (momentum, heat, freshwater and runoff)
289      !!
290      !! ** Method  : (1) READ each fluxes in NetCDF files:
291      !!      the 10m wind velocity (i-component) (m/s)    at T-point
292      !!      the 10m wind velocity (j-component) (m/s)    at T-point
293      !!      the 10m or 2m specific humidity     ( % )
294      !!      the solar heat                      (W/m2)
295      !!      the Long wave                       (W/m2)
296      !!      the 10m or 2m air temperature       (Kelvin)
297      !!      the total precipitation (rain+snow) (Kg/m2/s)
298      !!      the snow (solid prcipitation)       (kg/m2/s)
299      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T)
300      !!              (2) CALL blk_oce
301      !!
302      !!      C A U T I O N : never mask the surface stress fields
303      !!                      the stress is assumed to be in the (i,j) mesh referential
304      !!
305      !! ** Action  :   defined at each time-step at the air-sea interface
306      !!              - utau, vtau  i- and j-component of the wind stress
307      !!              - taum        wind stress module at T-point
308      !!              - wndm        wind speed  module at T-point over free ocean or leads in presence of sea-ice
309      !!              - qns, qsr    non-solar and solar heat fluxes
310      !!              - emp         upward mass flux (evapo. - precip.)
311      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present)
312      !!                            (set in limsbc(_2).F90)
313      !!
314      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
315      !!                   Brodeau et al. Ocean Modelling 2010
316      !!----------------------------------------------------------------------
317      INTEGER, INTENT(in) ::   kt   ! ocean time step
318      !!---------------------------------------------------------------------
319      !
320      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
321      !
322      !                                            ! compute the surface ocean fluxes using bulk formulea
323      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
324
325#if defined key_cice
326      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
327         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
328         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
329         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
330         ENDIF
331         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
332         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
333         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
334         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
335         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
336         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
337      ENDIF
338#endif
339      !
340   END SUBROUTINE sbc_blk
341
342
343   SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
344      !!---------------------------------------------------------------------
345      !!                     ***  ROUTINE blk_oce  ***
346      !!
347      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
348      !!      the ocean surface at each time step
349      !!
350      !! ** Method  :   bulk formulea for the ocean using atmospheric
351      !!      fields read in sbc_read
352      !!
353      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
354      !!              - vtau    : j-component of the stress at V-point  (N/m2)
355      !!              - taum    : Wind stress module at T-point         (N/m2)
356      !!              - wndm    : Wind speed module at T-point          (m/s)
357      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
358      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
359      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
360      !!
361      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
362      !!---------------------------------------------------------------------
363      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
364      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
365      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
366      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
367      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
368      !
369      INTEGER  ::   ji, jj               ! dummy loop indices
370      REAL(wp) ::   zztmp                ! local variable
371      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point
372      REAL(wp), DIMENSION(:,:), POINTER ::   zsq               ! specific humidity at pst
373      REAL(wp), DIMENSION(:,:), POINTER ::   zqlw, zqsb        ! long wave and sensible heat fluxes
374      REAL(wp), DIMENSION(:,:), POINTER ::   zqla, zevap       ! latent heat fluxes and evaporation
375      REAL(wp), DIMENSION(:,:), POINTER ::   zst               ! surface temperature in Kelvin
376      REAL(wp), DIMENSION(:,:), POINTER ::   zU_zu             ! bulk wind speed at height zu  [m/s]
377      REAL(wp), DIMENSION(:,:), POINTER ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
378      REAL(wp), DIMENSION(:,:), POINTER ::   zrhoa             ! density of air   [kg/m^3]
379      !!---------------------------------------------------------------------
380      !
381      IF( nn_timing == 1 )  CALL timing_start('blk_oce')
382      !
383      CALL wrk_alloc( jpi,jpj,   zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap )
384      CALL wrk_alloc( jpi,jpj,   zst, zU_zu, ztpot, zrhoa )
385      !
386
387      ! local scalars ( place there for vector optimisation purposes)
388      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
389
390      ! ----------------------------------------------------------------------------- !
391      !      0   Wind components and module at T-point relative to the moving ocean   !
392      ! ----------------------------------------------------------------------------- !
393
394      ! ... components ( U10m - U_oce ) at T-point (unmasked)
395!!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
396      zwnd_i(:,:) = 0._wp
397      zwnd_j(:,:) = 0._wp
398#if defined key_cyclone
399      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
400      DO jj = 2, jpjm1
401         DO ji = fs_2, fs_jpim1   ! vect. opt.
402            sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
403            sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
404         END DO
405      END DO
406#endif
407      DO jj = 2, jpjm1
408         DO ji = fs_2, fs_jpim1   ! vect. opt.
409            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
410            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
411         END DO
412      END DO
413      CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
414      CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
415      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
416      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
417         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
418
419      ! ----------------------------------------------------------------------------- !
420      !      I   Radiative FLUXES                                                     !
421      ! ----------------------------------------------------------------------------- !
422
423      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
424      zztmp = 1. - albo
425      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
426      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
427      ENDIF
428
429      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
430
431      ! ----------------------------------------------------------------------------- !
432      !     II    Turbulent FLUXES                                                    !
433      ! ----------------------------------------------------------------------------- !
434
435      ! ... specific humidity at SST and IST tmask(
436      zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
437      !!
438      !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
439      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
440      !!    (since reanalysis products provide T at z, not theta !)
441      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
442
443      SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
444      !
445      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2
446         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
447      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0
448         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
449      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
450         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
451      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF
452         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
453      CASE DEFAULT
454         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
455      END SELECT
456
457      !                          ! Compute true air density :
458      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
459         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
460      ELSE                                      ! At zt:
461         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
462      END IF
463
464!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
465!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
466
467      DO jj = 1, jpj             ! tau module, i and j component
468         DO ji = 1, jpi
469            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
470            taum  (ji,jj) = zztmp * wndm  (ji,jj)
471            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
472            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
473         END DO
474      END DO
475
476      !                          ! add the HF tau contribution to the wind stress module
477      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
478
479      CALL iom_put( "taum_oce", taum )   ! output wind stress module
480
481      ! ... utau, vtau at U- and V_points, resp.
482      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
483      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
484      DO jj = 1, jpjm1
485         DO ji = 1, fs_jpim1
486            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
487               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
488            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
489               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
490         END DO
491      END DO
492      CALL lbc_lnk( utau(:,:), 'U', -1. )
493      CALL lbc_lnk( vtau(:,:), 'V', -1. )
494
495
496      !  Turbulent fluxes over ocean
497      ! -----------------------------
498
499      ! zqla used as temporary array, for rho*U (common term of bulk formulae):
500      zqla(:,:) = zrhoa(:,:) * zU_zu(:,:)
501
502      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
503         !! q_air and t_air are given at 10m (wind reference height)
504         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
505         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
506      ELSE
507         !! q_air and t_air are not given at 10m (wind reference height)
508         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
509         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
510         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
511      ENDIF
512
513      zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
514
515
516      IF(ln_ctl) THEN
517         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
518         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
519         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
520         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
521         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
522            &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask )
523         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ')
524         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ')
525      ENDIF
526
527      ! ----------------------------------------------------------------------------- !
528      !     III    Total FLUXES                                                       !
529      ! ----------------------------------------------------------------------------- !
530      !
531      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
532         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
533      !
534      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
535         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
536         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
537         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
538         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
539         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
540         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1)
541      !
542#if defined key_lim3
543      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by LIM3)
544      qsr_oce(:,:) = qsr(:,:)
545#endif
546      !
547      IF ( nn_ice == 0 ) THEN
548         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
549         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
550         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
551         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
552         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
553         CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
554         CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
555         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac   ! output total precipitation [kg/m2/s]
556         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac   ! output solid precipitation [kg/m2/s]
557         CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow
558         CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation
559      ENDIF
560      !
561      IF(ln_ctl) THEN
562         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
563         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
564         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
565         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
566            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
567      ENDIF
568      !
569      CALL wrk_dealloc( jpi,jpj,   zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap )
570      CALL wrk_dealloc( jpi,jpj,   zst, zU_zu, ztpot, zrhoa )
571      !
572      IF( nn_timing == 1 )  CALL timing_stop('blk_oce')
573      !
574   END SUBROUTINE blk_oce
575
576#if defined key_lim3
577
578   SUBROUTINE blk_ice_tau
579      !!---------------------------------------------------------------------
580      !!                     ***  ROUTINE blk_ice_tau  ***
581      !!
582      !! ** Purpose :   provide the surface boundary condition over sea-ice
583      !!
584      !! ** Method  :   compute momentum using bulk formulation
585      !!                formulea, ice variables and read atmospheric fields.
586      !!                NB: ice drag coefficient is assumed to be a constant
587      !!---------------------------------------------------------------------
588      INTEGER  ::   ji, jj    ! dummy loop indices
589      !
590      REAL(wp), DIMENSION(:,:)  , POINTER :: zrhoa
591      !
592      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point
593      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point
594      !!---------------------------------------------------------------------
595      !
596      IF( nn_timing == 1 )  CALL timing_start('blk_ice_tau')
597      !
598      CALL wrk_alloc( jpi,jpj, zrhoa )
599
600      ! set transfer coefficients to default sea-ice values
601      Cd_atm(:,:) = Cd_ice
602      Ch_atm(:,:) = Cd_ice
603      Ce_atm(:,:) = Cd_ice
604
605      wndm_ice(:,:) = 0._wp      !!gm brutal....
606
607      ! ------------------------------------------------------------ !
608      !    Wind module relative to the moving ice ( U10m - U_ice )   !
609      ! ------------------------------------------------------------ !
610      SELECT CASE( cp_ice_msh )
611      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
612         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked)
613         DO jj = 2, jpjm1
614            DO ji = 2, jpim1   ! B grid : NO vector opt
615               ! ... scalar wind at T-point (fld being at T-point)
616               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   &
617                  &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  )
618               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   &
619                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  )
620               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
621            END DO
622         END DO
623         CALL lbc_lnk( wndm_ice, 'T',  1. )
624         !
625      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean)
626         DO jj = 2, jpjm1
627            DO ji = fs_2, fs_jpim1   ! vect. opt.
628               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
629               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
630               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
631            END DO
632         END DO
633         CALL lbc_lnk( wndm_ice, 'T',  1. )
634         !
635      END SELECT
636
637      ! Make ice-atm. drag dependent on ice concentration
638      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
639         CALL Cdn10_Lupkes2012( Cd_atm )
640         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
641      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
642         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 
643      ENDIF
644
645!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
646!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
647
648      ! local scalars ( place there for vector optimisation purposes)
649      ! Computing density of air! Way denser that 1.2 over sea-ice !!!
650      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
651
652      !!gm brutal....
653      utau_ice  (:,:) = 0._wp
654      vtau_ice  (:,:) = 0._wp
655      !!gm end
656
657      ! ------------------------------------------------------------ !
658      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
659      ! ------------------------------------------------------------ !
660      SELECT CASE( cp_ice_msh )
661      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
662         DO jj = 2, jpjm1
663            DO ji = 2, jpim1   ! B grid : NO vector opt
664               ! ... scalar wind at I-point (fld being at T-point)
665               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   &
666                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj)
667               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   &
668                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj)
669               ! ... ice stress at I-point
670               zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
671               utau_ice(ji,jj) = zwnorm_f * zwndi_f
672               vtau_ice(ji,jj) = zwnorm_f * zwndj_f
673            END DO
674         END DO
675         CALL lbc_lnk( utau_ice, 'I', -1. )
676         CALL lbc_lnk( vtau_ice, 'I', -1. )
677         !
678      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean)
679         DO jj = 2, jpjm1
680            DO ji = fs_2, fs_jpim1   ! vect. opt.
681               utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
682                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
683               vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
684                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
685            END DO
686         END DO
687         CALL lbc_lnk( utau_ice, 'U', -1. )
688         CALL lbc_lnk( vtau_ice, 'V', -1. )
689         !
690      END SELECT
691
692      IF(ln_ctl) THEN
693         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
694         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
695      ENDIF
696
697      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_tau')
698
699   END SUBROUTINE blk_ice_tau
700
701
702   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
703      !!---------------------------------------------------------------------
704      !!                     ***  ROUTINE blk_ice_flx  ***
705      !!
706      !! ** Purpose :   provide the surface boundary condition over sea-ice
707      !!
708      !! ** Method  :   compute heat and freshwater exchanged
709      !!                between atmosphere and sea-ice using bulk formulation
710      !!                formulea, ice variables and read atmmospheric fields.
711      !!
712      !! caution : the net upward water flux has with mm/day unit
713      !!---------------------------------------------------------------------
714      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
715      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
716      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
717      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
718      !!
719      INTEGER  ::   ji, jj, jl               ! dummy loop indices
720      REAL(wp) ::   zst2, zst3               ! local variable
721      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
722      REAL(wp) ::   zztmp, z1_lsub           !   -      -
723      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice
724      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables
725      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice
726      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice
727      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw        ! long wave heat sensitivity over ice
728      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb        ! sensible  heat sensitivity over ice
729      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3)
730      REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa
731
732      !!---------------------------------------------------------------------
733      !
734      IF( nn_timing == 1 )  CALL timing_start('blk_ice_flx')
735      !
736      CALL wrk_alloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb )
737      CALL wrk_alloc( jpi,jpj,       zrhoa )
738      !
739      ! local scalars
740      zcoef_dqlw   = 4.0 * 0.95 * Stef
741      zcoef_dqla   = -Ls * 11637800. * (-5897.8)
742      !
743      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
744      !
745      zztmp = 1. / ( 1. - albo )
746      !                                     ! ========================== !
747      DO jl = 1, jpl                        !  Loop over ice categories  !
748         !                                  ! ========================== !
749         DO jj = 1 , jpj
750            DO ji = 1, jpi
751               ! ----------------------------!
752               !      I   Radiative FLUXES   !
753               ! ----------------------------!
754               zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
755               zst3 = ptsu(ji,jj,jl) * zst2
756               ! Short Wave (sw)
757               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
758               ! Long  Wave (lw)
759               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
760               ! lw sensitivity
761               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
762
763               ! ----------------------------!
764               !     II    Turbulent FLUXES  !
765               ! ----------------------------!
766
767               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
768               ! Sensible Heat
769               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
770               ! Latent Heat
771               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
772                  &                ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
773               ! Latent heat sensitivity for ice (Dqla/Dt)
774               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
775                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) / zst2 * EXP(-5897.8 / ptsu(ji,jj,jl))
776               ELSE
777                  dqla_ice(ji,jj,jl) = 0._wp
778               ENDIF
779
780               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
781               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
782
783               ! ----------------------------!
784               !     III    Total FLUXES     !
785               ! ----------------------------!
786               ! Downward Non Solar flux
787               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
788               ! Total non solar heat flux sensitivity for ice
789               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
790            END DO
791            !
792         END DO
793         !
794      END DO
795      !
796      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s]
797      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s]
798      CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation
799      CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation
800
801      CALL wrk_alloc( jpi,jpj,   zevap, zsnw )
802
803      ! --- evaporation --- !
804      z1_lsub = 1._wp / Lsub
805      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation
806      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT
807      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean
808
809      ! --- evaporation minus precipitation --- !
810      zsnw(:,:) = 0._wp
811      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
812      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
813      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
814      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
815
816      ! --- heat flux associated with emp --- !
817      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
818         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
819         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
820         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
821      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
822         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
823
824      ! --- total solar and non solar fluxes --- !
825      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
826         &           + qemp_ice(:,:) + qemp_oce(:,:)
827      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
828
829      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
830      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
831
832      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
833      DO jl = 1, jpl
834         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
835                                   ! But we do not have Tice => consider it at 0degC => evap=0
836      END DO
837
838      CALL wrk_dealloc( jpi,jpj,   zevap, zsnw )
839
840      ! --- absorbed and transmitted shortwave radiation (W/m2) --- !
841      !
842      ! former coding was
843      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
844      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
845
846      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- !
847      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value
848      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1
849
850      qsr_ice_tr(:,:,:) = 0._wp
851
852      DO jl = 1, jpl
853         DO jj = 1, jpj
854            DO ji = 1, jpi
855
856               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1
857               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
858
859               IF ( phs(ji,jj,jl) <= 0.0_wp ) THEN
860                  zfrqsr_tr_i  = zfr1            + zfac * zfr2
861               ELSE
862                  zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque
863               ENDIF
864
865               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation
866
867            END DO
868         END DO
869      END DO
870                                                                                     
871
872      IF(ln_ctl) THEN
873         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
874         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
875         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
876         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
877         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
878         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
879      ENDIF
880
881      CALL wrk_dealloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb )
882      CALL wrk_dealloc( jpi,jpj,       zrhoa )
883      !
884      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_flx')
885
886   END SUBROUTINE blk_ice_flx
887   
888   
889
890   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi )
891   
892      !!---------------------------------------------------------------------
893      !!                     ***  ROUTINE blk_ice_qcn  ***
894      !!
895      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
896      !!                to force sea ice / snow thermodynamics
897      !!                in the case JULES coupler is emulated
898      !!               
899      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
900      !!                following the 0-layer Semtner (1976) approach
901      !!
902      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
903      !!              - qcn_ice : surface inner conduction flux (W/m2)
904      !!
905      !!---------------------------------------------------------------------
906      !!
907      INTEGER, INTENT(in)                        ::   k_monocat                 ! single-category option
908
909      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   ptsu                      ! sea ice / snow surface temperature
910
911      REAL(wp), DIMENSION(:,:)  , INTENT(in)     ::   ptb                       ! sea ice base temperature
912      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phs                       ! snow thickness
913      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phi                       ! sea ice thickness
914     
915      INTEGER                                    ::   ji, jj, jl                ! dummy loop indices
916      INTEGER                                    ::   iter                      !
917      REAL(wp)                                   ::   zfac, zfac2, zfac3        ! dummy factors
918      REAL(wp)                                   ::   zkeff_h, ztsu             !
919      REAL(wp)                                   ::   zqc, zqnet                !
920      REAL(wp)                                   ::   zhe, zqa0                 !
921           
922      INTEGER , PARAMETER                        ::   nit = 10                  ! number of iterations
923      REAL(wp), PARAMETER                        ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
924
925      REAL(wp), DIMENSION(:,:,:), POINTER        ::   zgfac                    ! enhanced conduction factor
926     
927      !!---------------------------------------------------------------------
928
929      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn')
930      !
931      CALL wrk_alloc( jpi,jpj,jpl,   zgfac )
932
933      ! -------------------------------------!
934      !      I   Enhanced conduction factor  !
935      ! -------------------------------------!
936      !
937      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3)
938      ! Fichefet and Morales Maqueda, JGR 1997
939      !
940      zgfac(:,:,:) = 1._wp
941     
942      SELECT CASE ( k_monocat ) 
943     
944         CASE ( 1 , 3 )
945     
946            zfac    =   1._wp /  ( rn_cnd_s + rcdic )
947            zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon
948            zfac3   =   2._wp / zepsilon
949     
950            DO jl = 1, jpl               
951               DO jj = 1 , jpj
952                  DO ji = 1, jpi
953                     !                                ! Effective thickness
954                     zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac
955                     !                                ! Enhanced conduction factor
956                     IF( zhe >=  zfac2 )  &
957                        zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) )
958                  END DO
959               END DO
960            END DO
961     
962      END SELECT
963     
964      ! -------------------------------------------------------------!
965      !      II   Surface temperature and conduction flux            !
966      ! -------------------------------------------------------------!
967
968      zfac   =   rcdic * rn_cnd_s 
969                                      ! ========================== !
970      DO jl = 1, jpl                  !  Loop over ice categories  !
971         !                            ! ========================== !
972         DO jj = 1 , jpj
973            DO ji = 1, jpi
974               !                      ! Effective conductivity of the snow-ice system divided by thickness
975               zkeff_h   =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) )
976                                      ! Store initial surface temperature
977               ztsu      =   ptsu(ji,jj,jl)
978                                      ! Net initial atmospheric heat flux
979               zqa0      =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl)
980               
981               DO iter = 1, nit       ! --- Iteration loop
982               
983                   !                  ! Conduction heat flux through snow-ice system (>0 downwards)
984                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) )
985                   !                  ! Surface energy budget
986                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc
987                   !                  ! Temperature update
988                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )
989               
990               END DO
991               
992               ptsu(ji,jj,jl)    = MIN( rt0, ztsu )
993               
994               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
995               
996            END DO
997         END DO
998         
999      END DO
1000     
1001      CALL wrk_dealloc( jpi,jpj,jpl,   zgfac )
1002     
1003      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn')
1004
1005   END SUBROUTINE blk_ice_qcn
1006   
1007#endif
1008
1009   FUNCTION rho_air( ptak, pqa, pslp )
1010      !!-------------------------------------------------------------------------------
1011      !!                           ***  FUNCTION rho_air  ***
1012      !!
1013      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
1014      !!
1015      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
1016      !!-------------------------------------------------------------------------------
1017      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
1018      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
1019      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
1020      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
1021      !!-------------------------------------------------------------------------------
1022      !
1023      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
1024      !
1025   END FUNCTION rho_air
1026
1027
1028   FUNCTION cp_air( pqa )
1029      !!-------------------------------------------------------------------------------
1030      !!                           ***  FUNCTION cp_air  ***
1031      !!
1032      !! ** Purpose : Compute specific heat (Cp) of moist air
1033      !!
1034      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
1035      !!-------------------------------------------------------------------------------
1036      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
1037      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
1038      !!-------------------------------------------------------------------------------
1039      !
1040      Cp_air = Cp_dry + Cp_vap * pqa
1041      !
1042   END FUNCTION cp_air
1043
1044
1045   FUNCTION q_sat( ptak, pslp )
1046      !!----------------------------------------------------------------------------------
1047      !!                           ***  FUNCTION q_sat  ***
1048      !!
1049      !! ** Purpose : Specific humidity at saturation in [kg/kg]
1050      !!              Based on accurate estimate of "e_sat"
1051      !!              aka saturation water vapor (Goff, 1957)
1052      !!
1053      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
1054      !!----------------------------------------------------------------------------------
1055      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
1056      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
1057      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
1058      !
1059      INTEGER  ::   ji, jj         ! dummy loop indices
1060      REAL(wp) ::   ze_sat, ztmp   ! local scalar
1061      !!----------------------------------------------------------------------------------
1062      !
1063      DO jj = 1, jpj
1064         DO ji = 1, jpi
1065            !
1066            ztmp = rt0 / ptak(ji,jj)
1067            !
1068            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
1069            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
1070               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
1071               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
1072               !
1073            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
1074            !
1075         END DO
1076      END DO
1077      !
1078   END FUNCTION q_sat
1079
1080
1081   FUNCTION gamma_moist( ptak, pqa )
1082      !!----------------------------------------------------------------------------------
1083      !!                           ***  FUNCTION gamma_moist  ***
1084      !!
1085      !! ** Purpose : Compute the moist adiabatic lapse-rate.
1086      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
1087      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
1088      !!
1089      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
1090      !!----------------------------------------------------------------------------------
1091      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
1092      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
1093      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
1094      !
1095      INTEGER  ::   ji, jj         ! dummy loop indices
1096      REAL(wp) :: zrv, ziRT        ! local scalar
1097      !!----------------------------------------------------------------------------------
1098      !
1099      DO jj = 1, jpj
1100         DO ji = 1, jpi
1101            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
1102            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
1103            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) )
1104         END DO
1105      END DO
1106      !
1107   END FUNCTION gamma_moist
1108
1109
1110   FUNCTION L_vap( psst )
1111      !!---------------------------------------------------------------------------------
1112      !!                           ***  FUNCTION L_vap  ***
1113      !!
1114      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
1115      !!
1116      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
1117      !!----------------------------------------------------------------------------------
1118      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
1119      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
1120      !!----------------------------------------------------------------------------------
1121      !
1122      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
1123      !
1124   END FUNCTION L_vap
1125
1126
1127#if defined key_lim3
1128
1129   SUBROUTINE Cdn10_Lupkes2012( Cd )
1130      !!----------------------------------------------------------------------
1131      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1132      !!
1133      !! ** Purpose :    Recompute the ice-atm drag at 10m height to make
1134      !!                 it dependent on edges at leads, melt ponds and flows.
1135      !!                 After some approximations, this can be resumed to a dependency
1136      !!                 on ice concentration.
1137      !!               
1138      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
1139      !!                 with the highest level of approximation: level4, eq.(59)
1140      !!                 The generic drag over a cell partly covered by ice can be re-written as follows:
1141      !!
1142      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
1143      !!
1144      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
1145      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1146      !!                    A is the concentration of ice minus melt ponds (if any)
1147      !!
1148      !!                 This new drag has a parabolic shape (as a function of A) starting at
1149      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1150      !!                 and going down to Cdi(say 1.4e-3) for A=1
1151      !!
1152      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
1153      !!                 => see Lupkes et al (2013)
1154      !!
1155      !! ** References : Lupkes et al. JGR 2012 (theory)
1156      !!                 Lupkes et al. GRL 2013 (application to GCM)
1157      !!
1158      !!----------------------------------------------------------------------
1159      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1160      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
1161      REAL(wp), PARAMETER ::   znu   = 1._wp
1162      REAL(wp), PARAMETER ::   zmu   = 1._wp
1163      REAL(wp), PARAMETER ::   zbeta = 1._wp
1164      REAL(wp)            ::   zcoef
1165      !!----------------------------------------------------------------------
1166      zcoef = znu + 1._wp / ( 10._wp * zbeta )
1167
1168      ! generic drag over a cell partly covered by ice
1169      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
1170      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
1171      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
1172
1173      ! ice-atm drag
1174      Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag
1175         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
1176     
1177   END SUBROUTINE Cdn10_Lupkes2012
1178
1179
1180   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1181      !!----------------------------------------------------------------------
1182      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1183      !!
1184      !! ** pUrpose :    1lternative turbulent transfert coefficients formulation
1185      !!                 between sea-ice and atmosphere with distinct momentum
1186      !!                 and heat coefficients depending on sea-ice concentration
1187      !!                 and atmospheric stability (no meltponds effect for now).
1188      !!               
1189      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
1190      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1191      !!                 it considers specific skin and form drags (Andreas et al. 2010)
1192      !!                 to compute neutral transfert coefficients for both heat and
1193      !!                 momemtum fluxes. Atmospheric stability effect on transfert
1194      !!                 coefficient is also taken into account following Louis (1979).
1195      !!
1196      !! ** References : Lupkes et al. JGR 2015 (theory)
1197      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
1198      !!
1199      !!----------------------------------------------------------------------
1200      !
1201      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1202      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
1203      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat
1204      !
1205      ! ECHAM6 constants
1206      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
1207      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
1208      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
1209      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
1210      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
1211      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
1212      REAL(wp), PARAMETER ::   zc2          = zc * zc
1213      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
1214      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
1215      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
1216      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
1217      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
1218      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
1219      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
1220      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
1221      !
1222      INTEGER  ::   ji, jj         ! dummy loop indices
1223      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
1224      REAL(wp) ::   zrib_o, zrib_i
1225      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1226      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
1227      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
1228      REAL(wp) ::   zCdn_form_tmp
1229      !!----------------------------------------------------------------------
1230
1231      ! Momentum Neutral Transfert Coefficients (should be a constant)
1232      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
1233      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
1234      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details)
1235      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1236
1237      ! Heat Neutral Transfert Coefficients
1238      zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52 (cf Lupkes email for details)
1239     
1240      ! Atmospheric and Surface Variables
1241      zst(:,:)     = sst_m(:,:) + rt0                                       ! convert SST from Celcius to Kelvin
1242      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)  , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
1243      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
1244      !
1245      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
1246         DO ji = fs_2, fs_jpim1
1247            ! Virtual potential temperature [K]
1248            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
1249            zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
1250            zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
1251           
1252            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1253            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
1254            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
1255           
1256            ! Momentum and Heat Neutral Transfert Coefficients
1257            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
1258            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
1259                       
1260            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1261            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1262            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details)
1263            IF( zrib_o <= 0._wp ) THEN
1264               zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) )  ! Eq. 10
1265               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
1266                  &             )**zgamma )**z1_gamma
1267            ELSE
1268               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
1269               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
1270            ENDIF
1271           
1272            IF( zrib_i <= 0._wp ) THEN
1273               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
1274               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
1275            ELSE
1276               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
1277               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
1278            ENDIF
1279           
1280            ! Momentum Transfert Coefficients (Eq. 38)
1281            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
1282               &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
1283           
1284            ! Heat Transfert Coefficients (Eq. 49)
1285            Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
1286               &        zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
1287            !
1288         END DO
1289      END DO
1290      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. )
1291      !
1292   END SUBROUTINE Cdn10_Lupkes2015
1293
1294#endif
1295
1296   !!======================================================================
1297END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.