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 NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90 @ 12099

Last change on this file since 12099 was 12099, checked in by laurent, 5 years ago

Tiny fix for AGRIF to compile...

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