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

Context Navigation

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

Last change on this file since 11284 was 11284, checked in by laurent, 5 years ago
LB: gooodbye "coare3p5", hello "coare3p6" + cleaner air humidity handling + support for Relative humidity.
Property svn:keywords set to `Id`
File size: 70.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: