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/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90 @ 14022

Last change on this file since 14022 was 14022, checked in by laurent, 3 years ago
Caught up with trunk rev 14021...
Property svn:keywords set to `Id`
File size: 79.8 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	!! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE)
21	!! 4.2 ! 2020-12 (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
22	!!----------------------------------------------------------------------
23
24	!!----------------------------------------------------------------------
25	!! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition
26	!! sbc_blk : bulk formulation as ocean surface boundary condition
27	!! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE)
28	!! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE)
29	!! sea-ice case only :
30	!! blk_ice_1 : provide the air-ice stress
31	!! blk_ice_2 : provide the heat and mass fluxes at air-ice interface
32	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
33	!!----------------------------------------------------------------------
34	USE oce ! ocean dynamics and tracers
35	USE dom_oce ! ocean space and time domain
36	USE phycst ! physical constants
37	USE fldread ! read input fields
38	USE sbc_oce ! Surface boundary condition: ocean fields
39	USE cyclone ! Cyclone 10m wind form trac of cyclone centres
40	USE sbcdcy ! surface boundary condition: diurnal cycle
41	USE sbcwave , ONLY : cdn_wave ! wave module
42	USE lib_fortran ! to use key_nosignedzero
43	!
44	#if defined key_si3
45	USE sbc_ice ! Surface boundary condition: ice fields #LB? ok to be in 'key_si3' ???
46	USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif, nn_qtrice
47	USE icevar ! for CALL ice_var_snwblow
48	USE sbcblk_algo_ice_an05
49	USE sbcblk_algo_ice_lu12
50	USE sbcblk_algo_ice_lg15
51	#endif
52	USE sbcblk_algo_ncar ! => turb_ncar : NCAR - (formerly known as CORE, Large & Yeager, 2009)
53	USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003)
54	USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013)
55	USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 45r1)
56	USE sbcblk_algo_andreas ! => turb_andreas : Andreas et al. 2015
57	!
58	USE iom ! I/O manager library
59	USE in_out_manager ! I/O manager
60	USE lib_mpp ! distribued memory computing library
61	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
62	USE prtctl ! Print control
63
64	USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
65
66	IMPLICIT NONE
67	PRIVATE
68
69	PUBLIC sbc_blk_init ! called in sbcmod
70	PUBLIC sbc_blk ! called in sbcmod
71	PUBLIC blk_oce_1 ! called in sbcabl
72	PUBLIC blk_oce_2 ! called in sbcabl
73	#if defined key_si3
74	PUBLIC blk_ice_1 ! routine called in icesbc
75	PUBLIC blk_ice_2 ! routine called in icesbc
76	PUBLIC blk_ice_qcn ! routine called in icesbc
77	#endif
78
79	INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
80	INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
81	INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin)
82	INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % )
83	INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2)
84	INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2)
85	INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
86	INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s)
87	INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa)
88	INTEGER , PUBLIC, PARAMETER :: jp_uoatm = 10 ! index of surface current (i-component)
89	! ! seen by the atmospheric forcing (m/s) at T-point
90	INTEGER , PUBLIC, PARAMETER :: jp_voatm = 11 ! index of surface current (j-component)
91	! ! seen by the atmospheric forcing (m/s) at T-point
92	INTEGER , PUBLIC, PARAMETER :: jp_cc = 12 ! index of cloud cover (-) range:0-1
93	INTEGER , PUBLIC, PARAMETER :: jp_hpgi = 13 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
94	INTEGER , PUBLIC, PARAMETER :: jp_hpgj = 14 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
95	INTEGER , PUBLIC, PARAMETER :: jpfld = 14 ! maximum number of files to read
96
97	! Warning: keep this structure allocatable for Agrif...
98	TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read)
99
100	! !!* Namelist namsbc_blk : bulk parameters
101	LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008)
102	LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003)
103	LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013)
104	LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1)
105	LOGICAL :: ln_ANDREAS ! "ANDREAS" algorithm (Andreas et al. 2015)
106	!
107	!#LB:
108	LOGICAL :: ln_Cx_ice_cst ! use constant air-ice bulk transfer coefficients (value given in namelist's rn_Cd_i, rn_Ce_i & rn_Ch_i)
109	REAL(wp) :: rn_Cd_i, rn_Ce_i, rn_Ch_i ! values for " "
110	LOGICAL :: ln_Cx_ice_AN05 ! air-ice bulk transfer coefficients based on Andreas et al., 2005
111	LOGICAL :: ln_Cx_ice_LU12 ! air-ice bulk transfer coefficients based on Lupkes et al., 2012
112	LOGICAL :: ln_Cx_ice_LG15 ! air-ice bulk transfer coefficients based on Lupkes & Gryanik, 2015
113	!#LB.
114	!
115	LOGICAL :: ln_crt_fbk ! Add surface current feedback to the wind stress computation (Renault et al. 2020)
116	REAL(wp) :: rn_stau_a ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta
117	REAL(wp) :: rn_stau_b !
118	!
119	REAL(wp) :: rn_pfac ! multiplication factor for precipitation
120	REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation
121	REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
122	REAL(wp) :: rn_zu ! z(u) : height of wind measurements
123	!
124	INTEGER :: nn_iter_algo ! Number of iterations in bulk param. algo ("stable ABL + weak wind" requires more)
125
126	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: theta_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
127
128	#if defined key_si3
129	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice !#LB transfert coefficients over ice
130	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: theta_zu_i, q_zu_i !#LB fixme ! air temp. and spec. hum. over ice at wind speed height (L15 bulk scheme)
131	#endif
132
133
134	LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB
135	LOGICAL :: ln_skin_wl ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
136	LOGICAL :: ln_humi_sph ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB
137	LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB
138	LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB
139	LOGICAL :: ln_tair_pot ! temperature read in files ("sn_tair") is already potential temperature (not absolute)
140	!
141	INTEGER :: nhumi ! choice of the bulk algorithm
142	! ! associated indices:
143	INTEGER, PARAMETER :: np_humi_sph = 1
144	INTEGER, PARAMETER :: np_humi_dpt = 2
145	INTEGER, PARAMETER :: np_humi_rlh = 3
146
147	INTEGER :: nblk ! choice of the bulk algorithm
148	! ! associated indices:
149	INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008)
150	INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003)
151	INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013)
152	INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 45r1)
153	INTEGER, PARAMETER :: np_ANDREAS = 5 ! "ANDREAS" algorithm (Andreas et al. 2015)
154
155	!#LB:
156	#if defined key_si3
157	! Same, over sea-ice:
158	INTEGER :: nblk_ice ! choice of the bulk algorithm
159	! ! associated indices:
160	INTEGER, PARAMETER :: np_ice_cst = 1 ! constant transfer coefficients
161	INTEGER, PARAMETER :: np_ice_an05 = 2 ! Andreas et al., 2005
162	INTEGER, PARAMETER :: np_ice_lu12 = 3 ! Lupkes el al., 2012
163	INTEGER, PARAMETER :: np_ice_lg15 = 4 ! Lupkes & Gryanik, 2015
164	#endif
165	!LB.
166
167
168
169	!! * Substitutions
170	# include "do_loop_substitute.h90"
171	!!----------------------------------------------------------------------
172	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
173	!! $Id$
174	!! Software governed by the CeCILL license (see ./LICENSE)
175	!!----------------------------------------------------------------------
176	CONTAINS
177
178	INTEGER FUNCTION sbc_blk_alloc()
179	!!-------------------------------------------------------------------
180	!! * ROUTINE sbc_blk_alloc *
181	!!-------------------------------------------------------------------
182	ALLOCATE( theta_zu(jpi,jpj), q_zu(jpi,jpj), STAT=sbc_blk_alloc )
183	CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
184	IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
185	END FUNCTION sbc_blk_alloc
186
187	#if defined key_si3
188	INTEGER FUNCTION sbc_blk_ice_alloc()
189	!!-------------------------------------------------------------------
190	!! * ROUTINE sbc_blk_ice_alloc *
191	!!-------------------------------------------------------------------
192	ALLOCATE( Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), &
193	& theta_zu_i(jpi,jpj), q_zu_i(jpi,jpj), STAT=sbc_blk_ice_alloc )
194	CALL mpp_sum ( 'sbcblk', sbc_blk_ice_alloc )
195	IF( sbc_blk_ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_ice_alloc: failed to allocate arrays' )
196	END FUNCTION sbc_blk_ice_alloc
197	#endif
198
199
200	SUBROUTINE sbc_blk_init
201	!!---------------------------------------------------------------------
202	!! * ROUTINE sbc_blk_init *
203	!!
204	!! ** Purpose : choose and initialize a bulk formulae formulation
205	!!
206	!! ** Method :
207	!!
208	!!----------------------------------------------------------------------
209	INTEGER :: jfpr ! dummy loop indice and argument
210	INTEGER :: ios, ierror, ioptio ! Local integer
211	!!
212	CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files
213	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
214	TYPE(FLD_N) :: sn_wndi, sn_wndj , sn_humi, sn_qsr ! informations about the fields to be read
215	TYPE(FLD_N) :: sn_qlw , sn_tair , sn_prec, sn_snow ! " "
216	TYPE(FLD_N) :: sn_slp , sn_uoatm, sn_voatm ! " "
217	TYPE(FLD_N) :: sn_cc, sn_hpgi, sn_hpgj ! " "
218	INTEGER :: ipka ! number of levels in the atmospheric variable
219	NAMELIST/namsbc_blk/ ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, ln_ANDREAS, & ! bulk algorithm
220	& rn_zqt, rn_zu, nn_iter_algo, ln_skin_cs, ln_skin_wl, &
221	& rn_pfac, rn_efac, &
222	& ln_crt_fbk, rn_stau_a, rn_stau_b, & ! current feedback
223	& ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tair_pot, &
224	& ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i, &
225	& ln_Cx_ice_AN05, ln_Cx_ice_LU12, ln_Cx_ice_LG15, &
226	& cn_dir, &
227	& sn_wndi, sn_wndj, sn_qsr, sn_qlw , & ! input fields
228	& sn_tair, sn_humi, sn_prec, sn_snow, sn_slp, &
229	& sn_uoatm, sn_voatm, sn_cc, sn_hpgi, sn_hpgj
230
231	! cool-skin / warm-layer !LB
232	!!---------------------------------------------------------------------
233	!
234	! ! allocate sbc_blk_core array
235	IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
236	!
237	#if defined key_si3
238	IF( sbc_blk_ice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard ice arrays' )
239	#endif
240	!
241	! !** read bulk namelist
242	READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
243	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' )
244	!
245	READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
246	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' )
247	!
248	IF(lwm) WRITE( numond, namsbc_blk )
249	!
250	! !** initialization of the chosen bulk formulae (+ check)
251	! !* select the bulk chosen in the namelist and check the choice
252	ioptio = 0
253	IF( ln_NCAR ) THEN
254	nblk = np_NCAR ; ioptio = ioptio + 1
255	ENDIF
256	IF( ln_COARE_3p0 ) THEN
257	nblk = np_COARE_3p0 ; ioptio = ioptio + 1
258	ENDIF
259	IF( ln_COARE_3p6 ) THEN
260	nblk = np_COARE_3p6 ; ioptio = ioptio + 1
261	ENDIF
262	IF( ln_ECMWF ) THEN
263	nblk = np_ECMWF ; ioptio = ioptio + 1
264	ENDIF
265	IF( ln_ANDREAS ) THEN
266	nblk = np_ANDREAS ; ioptio = ioptio + 1
267	ENDIF
268	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
269
270	! !** initialization of the cool-skin / warm-layer parametrization
271	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
272	!! Some namelist sanity tests:
273	IF( ln_NCAR ) &
274	& CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' )
275	IF( ln_ANDREAS ) &
276	& CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with ANDREAS algorithm' )
277	IF( nn_fsbc /= 1 ) &
278	& CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
279	END IF
280
281	IF( ln_skin_wl ) THEN
282	!! Check if the frequency of downwelling solar flux input makes sense and if ln_dm2dc=T if it is daily!
283	IF( (sn_qsr%freqh < 0.).OR.(sn_qsr%freqh > 24.) ) &
284	& CALL ctl_stop( 'sbc_blk_init: Warm-layer param. (ln_skin_wl) not compatible with freq. of solar flux > daily' )
285	IF( (sn_qsr%freqh == 24.).AND.(.NOT. ln_dm2dc) ) &
286	& CALL ctl_stop( 'sbc_blk_init: Please set ln_dm2dc=T for warm-layer param. (ln_skin_wl) to work properly' )
287	END IF
288
289	ioptio = 0
290	IF( ln_humi_sph ) THEN
291	nhumi = np_humi_sph ; ioptio = ioptio + 1
292	ENDIF
293	IF( ln_humi_dpt ) THEN
294	nhumi = np_humi_dpt ; ioptio = ioptio + 1
295	ENDIF
296	IF( ln_humi_rlh ) THEN
297	nhumi = np_humi_rlh ; ioptio = ioptio + 1
298	ENDIF
299	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' )
300	!
301	IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr
302	IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
303	IF( sn_qsr%ln_tint ) THEN
304	CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', &
305	& ' ==> We force time interpolation = .false. for qsr' )
306	sn_qsr%ln_tint = .false.
307	ENDIF
308	ENDIF
309
310	#if defined key_si3
311	ioptio = 0
312	IF( ln_Cx_ice_cst ) THEN
313	nblk_ice = np_ice_cst ; ioptio = ioptio + 1
314	ENDIF
315	IF( ln_Cx_ice_AN05 ) THEN
316	nblk_ice = np_ice_an05 ; ioptio = ioptio + 1
317	ENDIF
318	IF( ln_Cx_ice_LU12 ) THEN
319	nblk_ice = np_ice_lu12 ; ioptio = ioptio + 1
320	ENDIF
321	IF( ln_Cx_ice_LG15 ) THEN
322	nblk_ice = np_ice_lg15 ; ioptio = ioptio + 1
323	ENDIF
324	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one ice-atm bulk algorithm' )
325	#endif
326
327
328	! !* set the bulk structure
329	! !- store namelist information in an array
330	!
331	slf_i(jp_wndi ) = sn_wndi ; slf_i(jp_wndj ) = sn_wndj
332	slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw
333	slf_i(jp_tair ) = sn_tair ; slf_i(jp_humi ) = sn_humi
334	slf_i(jp_prec ) = sn_prec ; slf_i(jp_snow ) = sn_snow
335	slf_i(jp_slp ) = sn_slp ; slf_i(jp_cc ) = sn_cc
336	slf_i(jp_uoatm) = sn_uoatm ; slf_i(jp_voatm) = sn_voatm
337	slf_i(jp_hpgi ) = sn_hpgi ; slf_i(jp_hpgj ) = sn_hpgj
338	!
339	IF( .NOT. ln_abl ) THEN ! force to not use jp_hpgi and jp_hpgj, should already be done in namelist_* but we never know...
340	slf_i(jp_hpgi)%clname = 'NOT USED'
341	slf_i(jp_hpgj)%clname = 'NOT USED'
342	ENDIF
343	!
344	! !- allocate the bulk structure
345	ALLOCATE( sf(jpfld), STAT=ierror )
346	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
347	!
348	! !- fill the bulk structure with namelist informations
349	CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
350	!
351	DO jfpr= 1, jpfld
352	!
353	IF( ln_abl .AND. &
354	& ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR. &
355	& jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair ) ) THEN
356	ipka = jpka ! ABL: some fields are 3D input
357	ELSE
358	ipka = 1
359	ENDIF
360	!
361	ALLOCATE( sf(jfpr)%fnow(jpi,jpj,ipka) )
362	!
363	IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to default)
364	IF( jfpr == jp_slp ) THEN
365	sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp ! use standard pressure in Pa
366	ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN
367	sf(jfpr)%fnow(:,:,1:ipka) = 0._wp ! no precip or no snow or no surface currents
368	ELSEIF( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) THEN
369	IF( .NOT. ln_abl ) THEN
370	DEALLOCATE( sf(jfpr)%fnow ) ! deallocate as not used in this case
371	ELSE
372	sf(jfpr)%fnow(:,:,1:ipka) = 0._wp
373	ENDIF
374	ELSEIF( jfpr == jp_cc ) THEN
375	sf(jp_cc)%fnow(:,:,1:ipka) = pp_cldf
376	ELSE
377	WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr
378	CALL ctl_stop( ctmp1 )
379	ENDIF
380	ELSE !-- used field --!
381	IF( sf(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) ) ! allocate array for temporal interpolation
382	!
383	IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 ) &
384	& CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', &
385	& ' This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' )
386	ENDIF
387	END DO
388	!
389	IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient
390	rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level
391	rn_zu = ght_abl(2)
392	IF(lwp) WRITE(numout,*)
393	IF(lwp) WRITE(numout,*) ' ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude'
394	ENDIF
395	!
396	!
397	IF(lwp) THEN !** Control print
398	!
399	WRITE(numout,) ! namelist
400	WRITE(numout,*) ' Namelist namsbc_blk (other than data information):'
401	WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR
402	WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0
403	WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013) ln_COARE_3p6 = ', ln_COARE_3p6
404	WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF
405	WRITE(numout,*) ' "ANDREAS" algorithm (Andreas et al. 2015) ln_ANDREAS = ', ln_ANDREAS
406	WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt
407	WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
408	WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
409	WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
410	WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
411	WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk
412	IF(ln_crt_fbk) THEN
413	WRITE(numout,) ' Renault et al. 2020, eq. 10: Stau = Alpha Wnd + Beta'
414	WRITE(numout,*) ' Alpha rn_stau_a = ', rn_stau_a
415	WRITE(numout,*) ' Beta rn_stau_b = ', rn_stau_b
416	ENDIF
417	!
418	WRITE(numout,*)
419	SELECT CASE( nblk ) !* Print the choice of bulk algorithm
420	CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)'
421	CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)'
422	CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)'
423	CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)'
424	CASE( np_ANDREAS ) ; WRITE(numout,*) ' ==>>> "ANDREAS" algorithm (Andreas et al. 2015)'
425	END SELECT
426	!
427	WRITE(numout,*)
428	WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs
429	WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl
430	!
431	WRITE(numout,*)
432	SELECT CASE( nhumi ) !* Print the choice of air humidity
433	CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]'
434	CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]'
435	CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]'
436	END SELECT
437	!
438	!#LB:
439	#if defined key_si3
440	IF( nn_ice > 0 ) THEN
441	WRITE(numout,*)
442	WRITE(numout,*) ' use constant ice-atm bulk transfer coeff. ln_Cx_ice_cst = ', ln_Cx_ice_cst
443	WRITE(numout,*) ' use ice-atm bulk coeff. from Andreas et al., 2005 ln_Cx_ice_AN05 = ', ln_Cx_ice_AN05
444	WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes et al., 2012 ln_Cx_ice_LU12 = ', ln_Cx_ice_LU12
445	WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes & Gryanik, 2015 ln_Cx_ice_LG15 = ', ln_Cx_ice_LG15
446	ENDIF
447	WRITE(numout,*)
448	SELECT CASE( nblk_ice ) !* Print the choice of bulk algorithm
449	CASE( np_ice_cst )
450	WRITE(numout,*) ' ==>>> Constant bulk transfer coefficients over sea-ice:'
451	WRITE(numout,*) ' => Cd_ice, Ce_ice, Ch_ice =', REAL(rn_Cd_i,4), REAL(rn_Ce_i,4), REAL(rn_Ch_i,4)
452	IF( (rn_Cd_i<0._wp).OR.(rn_Cd_i>1.E-2_wp).OR.(rn_Ce_i<0._wp).OR.(rn_Ce_i>1.E-2_wp).OR.(rn_Ch_i<0._wp).OR.(rn_Ch_i>1.E-2_wp) ) &
453	& CALL ctl_stop( 'Be realistic in your pick of Cd_ice, Ce_ice & Ch_ice ! (0 < Cx < 1.E-2)')
454	CASE( np_ice_an05 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Andreas et al, 2005'
455	CASE( np_ice_lu12 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes et al, 2012'
456	CASE( np_ice_lg15 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes & Gryanik, 2015'
457	END SELECT
458	#endif
459	!#LB.
460	!
461	ENDIF
462	!
463	END SUBROUTINE sbc_blk_init
464
465
466	SUBROUTINE sbc_blk( kt )
467	!!---------------------------------------------------------------------
468	!! * ROUTINE sbc_blk *
469	!!
470	!! ** Purpose : provide at each time step the surface ocean fluxes
471	!! (momentum, heat, freshwater and runoff)
472	!!
473	!! ** Method :
474	!! (1) READ each fluxes in NetCDF files:
475	!! the wind velocity (i-component) at z=rn_zu (m/s) at T-point
476	!! the wind velocity (j-component) at z=rn_zu (m/s) at T-point
477	!! the specific humidity at z=rn_zqt (kg/kg)
478	!! the air temperature at z=rn_zqt (Kelvin)
479	!! the solar heat (W/m2)
480	!! the Long wave (W/m2)
481	!! the total precipitation (rain+snow) (Kg/m2/s)
482	!! the snow (solid precipitation) (kg/m2/s)
483	!! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds)
484	!! (2) CALL blk_oce_1 and blk_oce_2
485	!!
486	!! C A U T I O N : never mask the surface stress fields
487	!! the stress is assumed to be in the (i,j) mesh referential
488	!!
489	!! ** Action : defined at each time-step at the air-sea interface
490	!! - utau, vtau i- and j-component of the wind stress
491	!! - taum wind stress module at T-point
492	!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
493	!! - qns, qsr non-solar and solar heat fluxes
494	!! - emp upward mass flux (evapo. - precip.)
495	!! - sfx salt flux due to freezing/melting (non-zero only if ice is present)
496	!!
497	!! ** References : Large & Yeager, 2004 / Large & Yeager, 2008
498	!! Brodeau et al. Ocean Modelling 2010
499	!!----------------------------------------------------------------------
500	INTEGER, INTENT(in) :: kt ! ocean time step
501	!!----------------------------------------------------------------------
502	REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp
503	REAL(wp) :: ztmp
504	!!----------------------------------------------------------------------
505	!
506	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
507
508	! Sanity/consistence test on humidity at first time step to detect potential screw-up:
509	IF( kt == nit000 ) THEN
510	IF(lwp) WRITE(numout,*) ''
511	#if defined key_agrif
512	IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ==='
513	#else
514	ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain
515	IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points!
516	ztmp = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztmp ! mean humidity over ocean on proc
517	SELECT CASE( nhumi )
518	CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!)
519	IF( (ztmp < 0._wp) .OR. (ztmp > 0.065) ) ztmp = -1._wp
520	CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K]
521	IF( (ztmp < 110._wp).OR.(ztmp > 320._wp) ) ztmp = -1._wp
522	CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%]
523	IF( (ztmp < 0._wp) .OR.(ztmp > 100._wp) ) ztmp = -1._wp
524	END SELECT
525	IF(ztmp < 0._wp) THEN
526	IF (lwp) WRITE(numout,'(" Mean humidity value found on proc #",i6.6," is: ",f10.5)') narea, ztmp
527	CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', &
528	& ' ==> check the unit in your input files' , &
529	& ' ==> check consistence of namelist choice: specific? relative? dew-point?', &
530	& ' ==> ln_humi_sph -> [kg/kg] \| ln_humi_rlh -> [%] \| ln_humi_dpt -> [K] !!!' )
531	END IF
532	END IF
533	IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ==='
534	#endif
535	IF(lwp) WRITE(numout,*) ''
536	END IF !IF( kt == nit000 )
537	! ! compute the surface ocean fluxes using bulk formulea
538	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
539
540	! Specific humidity of air at z=rn_zqt !
541	SELECT CASE( nhumi )
542	CASE( np_humi_sph )
543	q_air_zt(:,:) = sf(jp_humi )%fnow(:,:,1) ! what we read in file is already a spec. humidity!
544	CASE( np_humi_dpt )
545	IF((kt==nit000).AND.lwp) WRITE(numout,) ' ** sbc_blk() => computing q_air out of dew-point and P !'
546	q_air_zt(:,:) = q_sat( sf(jp_humi )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) )
547	CASE( np_humi_rlh )
548	IF((kt==nit000).AND.lwp) WRITE(numout,) ' ** sbc_blk() => computing q_air out of RH, t_air and slp !' !LBrm
549	q_air_zt(:,:) = q_air_rh( 0.01_wp*sf(jp_humi )%fnow(:,:,1), &
550	& sf(jp_tair )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) !#LB: 0.01 => RH is % percent in file
551	END SELECT
552
553	! POTENTIAL temperature of air at z=rn_zqt
554	! based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
555	! (most reanalysis products provide absolute temp., not potential temp.)
556	IF( ln_tair_pot ) THEN
557	! temperature read into file is already potential temperature, do nothing...
558	theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)
559	ELSE
560	! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...)
561	IF((kt==nit000).AND.lwp) WRITE(numout,) ' ** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!'
562	theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt
563	ENDIF
564	!
565	CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1), & ! <<= in
566	& theta_air_zt(:,:), q_air_zt(:,:), & ! <<= in
567	& sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in
568	& sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in
569	& sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs)
570	& tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out
571
572	CALL blk_oce_2( theta_air_zt(:,:), & ! <<= in
573	& sf(jp_qlw )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1), & ! <<= in
574	& sf(jp_snow )%fnow(:,:,1), tsk_m, & ! <<= in
575	& zsen, zlat, zevp ) ! <=> in out
576	ENDIF
577	!
578	#if defined key_cice
579	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
580	qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1)
581	IF( ln_dm2dc ) THEN
582	qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
583	ELSE
584	qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
585	ENDIF
586	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature instead ????
587	!tatm_ice(:,:) = theta_air_zt(:,:) !#LB: THIS! ?
588
589	qatm_ice(:,:) = q_air_zt(:,:) !#LB:
590
591	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
592	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
593	wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
594	wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
595	ENDIF
596	#endif
597	!
598	END SUBROUTINE sbc_blk
599
600
601	SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! inp
602	& pslp , pst , pu , pv, & ! inp
603	& puatm, pvatm, pdqsr , pdqlw , & ! inp
604	& ptsk , pssq , pcd_du, psen, plat, pevp ) ! out
605	!!---------------------------------------------------------------------
606	!! * ROUTINE blk_oce_1 *
607	!!
608	!! ** Purpose : if ln_blk=T, computes surface momentum, heat and freshwater fluxes
609	!! if ln_abl=T, computes Cd x \|U\|, Ch x \|U\|, Ce x \|U\| for ABL integration
610	!!
611	!! ** Method : bulk formulae using atmospheric fields from :
612	!! if ln_blk=T, atmospheric fields read in sbc_read
613	!! if ln_abl=T, the ABL model at previous time-step
614	!!
615	!! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg)
616	!! - pcd_du : Cd x \|dU\| at T-points (m/s)
617	!! - psen : sensible heat flux (W/m^2)
618	!! - plat : latent heat flux (W/m^2)
619	!! - pevp : evaporation (mm/s) #lolo
620	!!---------------------------------------------------------------------
621	INTEGER , INTENT(in ) :: kt ! time step index
622	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s]
623	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s]
624	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqair ! specific humidity at T-points [kg/kg]
625	REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin]
626	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa]
627	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celsius]
628	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s]
629	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s]
630	REAL(wp), INTENT(in ), DIMENSION(:,:) :: puatm ! surface current seen by the atm at T-point (i-component) [m/s]
631	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pvatm ! surface current seen by the atm at T-point (j-component) [m/s]
632	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqsr ! downwelling solar (shortwave) radiation at surface [W/m^2]
633	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2]
634	REAL(wp), INTENT( out), DIMENSION(:,:) :: ptsk ! skin temp. (or SST if CS & WL not used) [Celsius]
635	REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg]
636	REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du
637	REAL(wp), INTENT( out), DIMENSION(:,:) :: psen
638	REAL(wp), INTENT( out), DIMENSION(:,:) :: plat
639	REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp
640	!
641	INTEGER :: ji, jj ! dummy loop indices
642	REAL(wp) :: zztmp ! local variable
643	REAL(wp) :: zstmax, zstau
644	#if defined key_cyclone
645	REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
646	#endif
647	REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point
648	REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
649	REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean
650	REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean
651	REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean
652	REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2
653	!!---------------------------------------------------------------------
654	!
655	! local scalars ( place there for vector optimisation purposes)
656	! ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K)
657	ptsk(:,:) = pst(:,:) + rt0 ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!)
658
659	! --- cloud cover --- !
660	cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
661
662	! ----------------------------------------------------------------------------- !
663	! 0 Wind components and module at T-point relative to the moving ocean !
664	! ----------------------------------------------------------------------------- !
665
666	! ... components ( U10m - U_oce ) at T-point (unmasked)
667	#if defined key_cyclone
668	zwnd_i(:,:) = 0._wp
669	zwnd_j(:,:) = 0._wp
670	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
671	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
672	zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
673	zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
674	! ... scalar wind at T-point (not masked)
675	wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) )
676	END_2D
677	#else
678	! ... scalar wind module at T-point (not masked)
679	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
680	wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
681	END_2D
682	#endif
683	! ----------------------------------------------------------------------------- !
684	! I Solar FLUX !
685	! ----------------------------------------------------------------------------- !
686
687	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
688	zztmp = 1. - albo
689	IF( ln_dm2dc ) THEN
690	qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
691	ELSE
692	qsr(:,:) = zztmp * pdqsr(:,:) * tmask(:,:,1)
693	ENDIF
694
695
696	! ----------------------------------------------------------------------------- !
697	! II Turbulent FLUXES !
698	! ----------------------------------------------------------------------------- !
699
700	! specific humidity at SST
701	pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
702
703	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
704	!! Backup "bulk SST" and associated spec. hum.
705	zztmp1(:,:) = ptsk(:,:)
706	zztmp2(:,:) = pssq(:,:)
707	ENDIF
708
709	!! Time to call the user-selected bulk parameterization for
710	!! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more...
711	SELECT CASE( nblk )
712
713	CASE( np_NCAR )
714	CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
715	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
716	& nb_iter=nn_iter_algo )
717	!
718	CASE( np_COARE_3p0 )
719	CALL turb_coare3p0( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
720	& ln_skin_cs, ln_skin_wl, &
721	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
722	& nb_iter=nn_iter_algo, &
723	& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
724	!
725	CASE( np_COARE_3p6 )
726	CALL turb_coare3p6( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
727	& ln_skin_cs, ln_skin_wl, &
728	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
729	& nb_iter=nn_iter_algo, &
730	& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
731	!
732	CASE( np_ECMWF )
733	CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
734	& ln_skin_cs, ln_skin_wl, &
735	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
736	& nb_iter=nn_iter_algo, &
737	& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
738	!
739	CASE( np_ANDREAS )
740	CALL turb_andreas ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
741	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
742	& nb_iter=nn_iter_algo )
743	!
744	CASE DEFAULT
745	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' )
746	!
747	END SELECT
748
749	IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1))
750	IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1))
751	IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1))
752	!! LB: mainly here for debugging purpose:
753	IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ptair-rt0) * tmask(:,:,1)) ! potential temperature at z=zt
754	IF( iom_use('q_zt') ) CALL iom_put("q_zt", pqair * tmask(:,:,1)) ! specific humidity "
755	IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (theta_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu
756	IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity "
757	IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0
758	IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu
759
760	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
761	!! ptsk and pssq have been updated!!!
762	!!
763	!! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq:
764	WHERE ( fr_i(:,:) > 0.001_wp )
765	! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
766	ptsk(:,:) = zztmp1(:,:)
767	pssq(:,:) = zztmp2(:,:)
768	END WHERE
769	END IF
770
771	! Turbulent fluxes over ocean => BULK_FORMULA @ sbc_phy.F90
772	! -------------------------------------------------------------
773
774	IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp
775	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
776	zztmp = zU_zu(ji,jj)
777	wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod
778	pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
779	psen(ji,jj) = zztmp * zch_oce(ji,jj)
780	pevp(ji,jj) = zztmp * zce_oce(ji,jj)
781	rhoa(ji,jj) = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) )
782	END_2D
783	ELSE !== BLK formulation ==! turbulent fluxes computation
784	CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
785	& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
786	& wndm(:,:), zU_zu(:,:), pslp(:,:), &
787	& taum(:,:), psen(:,:), plat(:,:), &
788	& pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac )
789
790	psen(:,:) = psen(:,:) * tmask(:,:,1)
791	plat(:,:) = plat(:,:) * tmask(:,:,1)
792	taum(:,:) = taum(:,:) * tmask(:,:,1)
793	pevp(:,:) = pevp(:,:) * tmask(:,:,1)
794
795	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
796	IF( wndm(ji,jj) > 0._wp ) THEN
797	zztmp = taum(ji,jj) / wndm(ji,jj)
798	#if defined key_cyclone
799	ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
800	ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
801	#else
802	ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
803	ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
804	#endif
805	ELSE
806	ztau_i(ji,jj) = 0._wp
807	ztau_j(ji,jj) = 0._wp
808	ENDIF
809	END_2D
810
811	IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
812	zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp ) ! set the max value of Stau corresponding to a wind of 3 m/s (<0)
813	DO_2D( 0, 1, 0, 1 ) ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop
814	zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) ! stau (<0) must be smaller than zstmax
815	ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) )
816	ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
817	taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) )
818	END_2D
819	ENDIF
820
821	! ... utau, vtau at U- and V_points, resp.
822	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
823	! Note that coastal wind stress is not used in the code... so this extra care has no effect
824	DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T
825	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj ) ) &
826	& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
827	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji ,jj+1) ) &
828	& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
829	END_2D
830
831	IF( ln_crt_fbk ) THEN
832	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. )
833	ELSE
834	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
835	ENDIF
836
837	! Saving open-ocean wind-stress (module and components) on T-points:
838	CALL iom_put( "taum_oce", taum(:,:)*tmask(:,:,1) ) ! output wind stress module
839	!#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau" (U-grid) and vtau" (V-grid) does the job in: [DYN/dynatf.F90])
840	CALL iom_put( "utau_oce", ztau_i(:,:)*tmask(:,:,1) ) ! utau at T-points!
841	CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) ) ! vtau at T-points!
842
843	IF(sn_cfctl%l_prtctl) THEN
844	CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ')
845	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ')
846	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, &
847	& tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
848	CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ')
849	ENDIF
850	!
851	ENDIF !IF( ln_abl )
852
853	ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius
854
855	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
856	CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius
857	CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference...
858	ENDIF
859	!
860	END SUBROUTINE blk_oce_1
861
862
863	SUBROUTINE blk_oce_2( ptair, pdqlw, pprec, psnow, & ! <<= in
864	& ptsk, psen, plat, pevp ) ! <<= in
865	!!---------------------------------------------------------------------
866	!! * ROUTINE blk_oce_2 *
867	!!
868	!! ** Purpose : finalize the momentum, heat and freshwater fluxes computation
869	!! at the ocean surface at each time step knowing Cd, Ch, Ce and
870	!! atmospheric variables (from ABL or external data)
871	!!
872	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
873	!! - vtau : j-component of the stress at V-point (N/m2)
874	!! - taum : Wind stress module at T-point (N/m2)
875	!! - wndm : Wind speed module at T-point (m/s)
876	!! - qsr : Solar heat flux over the ocean (W/m2)
877	!! - qns : Non Solar heat flux over the ocean (W/m2)
878	!! - emp : evaporation minus precipitation (kg/m2/s)
879	!!---------------------------------------------------------------------
880	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair ! potential temperature of air #LB: confirm!
881	REAL(wp), INTENT(in), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2]
882	REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec
883	REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow
884	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptsk ! SKIN surface temperature [Celsius]
885	REAL(wp), INTENT(in), DIMENSION(:,:) :: psen
886	REAL(wp), INTENT(in), DIMENSION(:,:) :: plat
887	REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp
888	!
889	INTEGER :: ji, jj ! dummy loop indices
890	REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable
891	REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! net long wave radiative heat flux
892	!!---------------------------------------------------------------------
893	!
894	! local scalars ( place there for vector optimisation purposes)
895
896	! ----------------------------------------------------------------------------- !
897	! III Net longwave radiative FLUX !
898	! ----------------------------------------------------------------------------- !
899	!! #LB: now moved after Turbulent fluxes because must use the skin temperature rather than bulk SST
900	!! (ptsk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
901	zqlw(:,:) = qlw_net( pdqlw(:,:), ptsk(:,:)+rt0 )
902
903	! ----------------------------------------------------------------------------- !
904	! IV Total FLUXES !
905	! ----------------------------------------------------------------------------- !
906	!
907	emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.)
908	& - pprec(:,:) * rn_pfac ) * tmask(:,:,1)
909	!
910	qns(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) & ! Downward Non Solar
911	& - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip
912	& - pevp(:,:) * ptsk(:,:) * rcp & ! remove evap heat content at SST
913	& + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair
914	& * ( ptair(:,:) - rt0 ) * rcp &
915	& + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
916	& * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi
917	qns(:,:) = qns(:,:) * tmask(:,:,1)
918	!
919	#if defined key_si3
920	qns_oce(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) ! non solar without emp (only needed by SI3)
921	qsr_oce(:,:) = qsr(:,:)
922	#endif
923	!
924	CALL iom_put( "rho_air" , rhoa*tmask(:,:,1) ) ! output air density [kg/m^3]
925	CALL iom_put( "evap_oce" , pevp ) ! evaporation
926	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
927	CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean
928	CALL iom_put( "qla_oce" , plat ) ! output downward latent heat over the ocean
929	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
930	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
931	CALL iom_put( 'snowpre', sprecip ) ! Snow
932	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
933	!
934	IF ( nn_ice == 0 ) THEN
935	CALL iom_put( "qemp_oce" , qns-zqlw-psen-plat ) ! output downward heat content of E-P over the ocean
936	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
937	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
938	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
939	ENDIF
940	!
941	IF(sn_cfctl%l_prtctl) THEN
942	CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ')
943	CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen : ' )
944	CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat : ' )
945	CALL prt_ctl(tab2d_1=qns , clinfo1=' blk_oce_2: qns : ' )
946	CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ')
947	ENDIF
948	!
949	END SUBROUTINE blk_oce_2
950
951
952	#if defined key_si3
953	!!----------------------------------------------------------------------
954	!! 'key_si3' SI3 sea-ice model
955	!!----------------------------------------------------------------------
956	!! blk_ice_1 : provide the air-ice stress
957	!! blk_ice_2 : provide the heat and mass fluxes at air-ice interface
958	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
959	!!----------------------------------------------------------------------
960
961	SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, pqair, pslp , puice, pvice, ptsui, & ! inputs
962	& putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs
963	!!---------------------------------------------------------------------
964	!! * ROUTINE blk_ice_1 *
965	!!
966	!! ** Purpose : provide the surface boundary condition over sea-ice
967	!!
968	!! ** Method : compute momentum using bulk formulation
969	!! formulea, ice variables and read atmospheric fields.
970	!! NB: ice drag coefficient is assumed to be a constant
971	!!---------------------------------------------------------------------
972	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa]
973	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s]
974	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s]
975	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s]
976	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric wind at T-point [m/s]
977	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s]
978	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! "
979	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K]
980	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk
981	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk
982	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl
983	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl
984	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl
985	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl
986	!
987	INTEGER :: ji, jj ! dummy loop indices
988	REAL(wp) :: zootm_su ! sea-ice surface mean temperature
989	REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
990	REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! temporary array
991	!!---------------------------------------------------------------------
992	!
993	! LB: ptsui is in K !!!
994	!
995	! ------------------------------------------------------------ !
996	! Wind module relative to the moving ice ( U10m - U_ice ) !
997	! ------------------------------------------------------------ !
998	! C-grid ice dynamics : U & V-points (same as ocean)
999	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1000	wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
1001	END_2D
1002	!
1003	! Make ice-atm. drag dependent on ice concentration
1004
1005
1006	SELECT CASE( nblk_ice )
1007
1008	CASE( np_ice_cst )
1009	! Constant bulk transfer coefficients over sea-ice:
1010	Cd_ice(:,:) = rn_Cd_i
1011	Ch_ice(:,:) = rn_Ch_i
1012	Ce_ice(:,:) = rn_Ce_i
1013	! no height adjustment, keeping zt values:
1014	theta_zu_i(:,:) = ptair(:,:)
1015	q_zu_i(:,:) = pqair(:,:)
1016
1017	CASE( np_ice_an05 ) ! calculate new drag from Lupkes(2015) equations
1018	ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
1019	CALL turb_ice_an05( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, &
1020	& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
1021	!!
1022	CASE( np_ice_lu12 )
1023	ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
1024	CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
1025	& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
1026	!!
1027	CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations
1028	ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
1029	CALL turb_ice_lg15( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
1030	& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
1031	!!
1032	END SELECT
1033
1034	IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ice').OR.iom_use('utau_ice').OR.iom_use('vtau_ice') ) &
1035	& ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice !
1036
1037	IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp)
1038	IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp)
1039	IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp)
1040
1041
1042	IF( ln_blk ) THEN
1043	! ---------------------------------------------------- !
1044	! Wind stress relative to nonmoving ice ( U10m ) !
1045	! ---------------------------------------------------- !
1046	! supress moving ice in wind stress computation as we don't know how to do it properly...
1047	DO_2D( 0, 1, 0, 1 ) ! at T point
1048	zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
1049	putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
1050	pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
1051	END_2D
1052
1053	!#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!!
1054	IF(iom_use('taum_ice')) CALL iom_put('taum_ice', SQRT( putauiputaui + pvtauipvtaui )*ztmp )
1055	!#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau_oi" (U-grid) and vtau_oi" (V-grid) does the job in: [ICE/icedyn_rhg_evp.F90])
1056	IF(iom_use('utau_ice')) CALL iom_put("utau_ice", putaui*ztmp) ! utau at T-points!
1057	IF(iom_use('vtau_ice')) CALL iom_put("vtau_ice", pvtaui*ztmp) ! vtau at T-points!
1058
1059	!
1060	DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean).
1061	!#LB: QUESTION?? so SI3 expects wind stress vector to be provided at U & V points? Not at T-points ?
1062	! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology
1063	zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) )
1064	zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) )
1065	putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) )
1066	pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) )
1067	END_2D
1068	CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
1069	!
1070	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' &
1071	& , tab2d_2=pvtaui , clinfo2=' pvtaui : ' )
1072	ELSE ! ln_abl
1073	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1074	pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj)
1075	pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
1076	pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
1077	END_2D
1078	!#LB:
1079	pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq !
1080	!#LB.
1081	ENDIF !IF( ln_blk )
1082	!
1083	IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
1084	!
1085	END SUBROUTINE blk_ice_1
1086
1087
1088	SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, pqair, pslp, pdqlw, pprec, psnow )
1089	!!---------------------------------------------------------------------
1090	!! * ROUTINE blk_ice_2 *
1091	!!
1092	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
1093	!!
1094	!! ** Method : compute heat and freshwater exchanged
1095	!! between atmosphere and sea-ice using bulk formulation
1096	!! formulea, ice variables and read atmmospheric fields.
1097	!!
1098	!! caution : the net upward water flux has with mm/day unit
1099	!!---------------------------------------------------------------------
1100	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K]
1101	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
1102	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
1103	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
1104	REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair ! potential temperature of air #LB: okay ???
1105	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqair ! specific humidity of air
1106	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp
1107	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pdqlw
1108	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec
1109	REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow
1110	!!
1111	INTEGER :: ji, jj, jl ! dummy loop indices
1112	REAL(wp) :: zst, zst3, zsq ! local variable
1113	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
1114	REAL(wp) :: zztmp, zzblk, zztmp1, z1_rLsub ! - -
1115	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
1116	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
1117	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
1118	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
1119	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
1120	REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2
1121	REAL(wp), DIMENSION(jpi,jpj) :: ztri
1122	!!---------------------------------------------------------------------
1123	!
1124	zcoef_dqlw = 4._wp * emiss_i * stefan ! local scalars
1125	!
1126
1127	zztmp = 1. / ( 1. - albo )
1128	dqla_ice(:,:,:) = 0._wp
1129
1130	! ! ========================== !
1131	DO jl = 1, jpl ! Loop over ice categories !
1132	! ! ========================== !
1133	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1134
1135	zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K]
1136	zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present
1137
1138	! ----------------------------!
1139	! I Radiative FLUXES !
1140	! ----------------------------!
1141	! Short Wave (sw)
1142	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
1143
1144	! Long Wave (lw)
1145	zst3 = zst * zst * zst
1146	z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1)
1147	! lw sensitivity
1148	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
1149
1150	! ----------------------------!
1151	! II Turbulent FLUXES !
1152	! ----------------------------!
1153
1154	! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
1155
1156	! Common term in bulk F. equations...
1157	zzblk = rhoa(ji,jj) * wndm_ice(ji,jj)
1158
1159	! Sensible Heat
1160	zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj)
1161	z_qsb (ji,jj,jl) = zztmp1 * (zst - theta_zu_i(ji,jj))
1162	z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice)
1163
1164	! Latent Heat
1165	zztmp1 = zzblk * rLsub * Ce_ice(ji,jj)
1166	qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ???
1167	IF(qla_ice(ji,jj,jl)>0._wp) dqla_ice(ji,jj,jl) = zztmp1*dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
1168	! !#LB: dq_sat_dt_ice() in "sbc_phy.F90"
1169	!#LB: without this unjustified "condensation sensure":
1170	!qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj))
1171	!dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
1172
1173
1174	! ----------------------------!
1175	! III Total FLUXES !
1176	! ----------------------------!
1177	! Downward Non Solar flux
1178	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
1179	! Total non solar heat flux sensitivity for ice
1180	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ????
1181
1182	END_2D
1183	!
1184	END DO
1185	!
1186	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]
1187	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]
1188	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
1189	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
1190
1191	! --- evaporation --- !
1192	z1_rLsub = 1._wp / rLsub
1193	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation
1194	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT
1195	zevap (:,:) = emp(:,:) + tprecip(:,:) ! evaporation over ocean !LB: removed rn_efac here, correct???
1196
1197	! --- evaporation minus precipitation --- !
1198	zsnw(:,:) = 0._wp
1199	CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
1200	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
1201	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
1202	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
1203
1204	! --- heat flux associated with emp --- !
1205	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
1206	& + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair
1207	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
1208	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1209	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
1210	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1211
1212	! --- total solar and non solar fluxes --- !
1213	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
1214	& + qemp_ice(:,:) + qemp_oce(:,:)
1215	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
1216
1217	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
1218	qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1219
1220	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
1221	DO jl = 1, jpl
1222	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) rcpi * tmask(:,:,1) )
1223	! ! But we do not have Tice => consider it at 0degC => evap=0
1224	END DO
1225
1226	! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- !
1227	IF( nn_qtrice == 0 ) THEN
1228	! formulation derived from Grenfell and Maykut (1977), where transmission rate
1229	! 1) depends on cloudiness
1230	! 2) is 0 when there is any snow
1231	! 3) tends to 1 for thin ice
1232	ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm
1233	DO jl = 1, jpl
1234	WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
1235	qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
1236	ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm
1237	qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
1238	ELSEWHERE ! zero when hs>0
1239	qtr_ice_top(:,:,jl) = 0._wp
1240	END WHERE
1241	ENDDO
1242	ELSEIF( nn_qtrice == 1 ) THEN
1243	! formulation is derived from the thesis of M. Lebrun (2019).
1244	! It represents the best fit using several sets of observations
1245	! It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
1246	qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:)
1247	ENDIF
1248	!
1249	IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
1250	ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )
1251	IF( iom_use('evap_ao_cea' ) ) CALL iom_put( 'evap_ao_cea' , ztmp(:,:) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
1252	IF( iom_use('hflx_evap_cea') ) CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) ) ! heat flux from evap (cell average)
1253	ENDIF
1254	IF( iom_use('hflx_rain_cea') ) THEN
1255	ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
1256	IF( iom_use('hflx_rain_cea') ) CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) ) ! heat flux from rain (cell average)
1257	ENDIF
1258	IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea') ) THEN
1259	WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 )
1260	ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
1261	ELSEWHERE
1262	ztmp(:,:) = rcp * sst_m(:,:)
1263	ENDWHERE
1264	ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )
1265	IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea' , ztmp2(:,:) ) ! heat flux from snow (cell average)
1266	IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
1267	IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) * zsnw(:,:) ) ! heat flux from snow (over ice)
1268	ENDIF
1269	!
1270	IF(sn_cfctl%l_prtctl) THEN
1271	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
1272	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
1273	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
1274	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
1275	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
1276	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
1277	ENDIF
1278
1279	!#LB:
1280	! air-ice heat flux components that are not written from ice_stp()@icestp.F90:
1281	IF( iom_use('qla_ice') ) CALL iom_put( 'qla_ice', SUM( - qla_ice * a_i_b, dim=3 ) ) !#LB: sign consistent with what's done for ocean
1282	IF( iom_use('qsb_ice') ) CALL iom_put( 'qsb_ice', SUM( - z_qsb * a_i_b, dim=3 ) ) !#LB: ==> negative => loss of heat for sea-ice
1283	IF( iom_use('qlw_ice') ) CALL iom_put( 'qlw_ice', SUM( z_qlw * a_i_b, dim=3 ) )
1284	!#LB.
1285
1286	END SUBROUTINE blk_ice_2
1287
1288
1289	SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
1290	!!---------------------------------------------------------------------
1291	!! * ROUTINE blk_ice_qcn *
1292	!!
1293	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
1294	!! to force sea ice / snow thermodynamics
1295	!! in the case conduction flux is emulated
1296	!!
1297	!! ** Method : compute surface energy balance assuming neglecting heat storage
1298	!! following the 0-layer Semtner (1976) approach
1299	!!
1300	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
1301	!! - qcn_ice : surface inner conduction flux (W/m2)
1302	!!
1303	!!---------------------------------------------------------------------
1304	LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option
1305	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
1306	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
1307	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
1308	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
1309	!
1310	INTEGER , PARAMETER :: nit = 10 ! number of iterations
1311	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
1312	!
1313	INTEGER :: ji, jj, jl ! dummy loop indices
1314	INTEGER :: iter ! local integer
1315	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
1316	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
1317	REAL(wp) :: zqc, zqnet !
1318	REAL(wp) :: zhe, zqa0 !
1319	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
1320	!!---------------------------------------------------------------------
1321
1322	! -------------------------------------!
1323	! I Enhanced conduction factor !
1324	! -------------------------------------!
1325	! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
1326	! Fichefet and Morales Maqueda, JGR 1997
1327	!
1328	zgfac(:,:,:) = 1._wp
1329
1330	IF( ld_virtual_itd ) THEN
1331	!
1332	zfac = 1._wp / ( rn_cnd_s + rcnd_i )
1333	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1334	zfac3 = 2._wp / zepsilon
1335	!
1336	DO jl = 1, jpl
1337	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1338	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness
1339	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
1340	END_2D
1341	END DO
1342	!
1343	ENDIF
1344
1345	! -------------------------------------------------------------!
1346	! II Surface temperature and conduction flux !
1347	! -------------------------------------------------------------!
1348	!
1349	zfac = rcnd_i * rn_cnd_s
1350	!
1351	DO jl = 1, jpl
1352	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1353	!
1354	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
1355	& ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1356	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
1357	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
1358	zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1359	!
1360	DO iter = 1, nit ! --- Iterative loop
1361	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
1362	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
1363	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
1364	END DO
1365	!
1366	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
1367	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1368	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1369	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) ) &
1370	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1371
1372	! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1373	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)
1374
1375	END_2D
1376	!
1377	END DO
1378	!
1379	END SUBROUTINE blk_ice_qcn
1380
1381	#endif
1382
1383	!!======================================================================
1384	END MODULE sbcblk

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

Download in other formats: