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 @ 13893

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

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

Download in other formats: