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

Last change on this file since 14021 was 14021, checked in by laurent, 3 years ago
Caught up with trunk rev 14020...
Property svn:keywords set to `Id`
File size: 79.6 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
59	!
60	USE iom ! I/O manager library
61	USE in_out_manager ! I/O manager
62	USE lib_mpp ! distribued memory computing library
63	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
64	USE prtctl ! Print control
65
66	USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
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_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient
392	rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level
393	rn_zu = ght_abl(2)
394	IF(lwp) WRITE(numout,*)
395	IF(lwp) WRITE(numout,*) ' ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude'
396	ENDIF
397	!
398	!
399	IF(lwp) THEN !** Control print
400	!
401	WRITE(numout,) ! namelist
402	WRITE(numout,*) ' Namelist namsbc_blk (other than data information):'
403	WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR
404	WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0
405	WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013) ln_COARE_3p6 = ', ln_COARE_3p6
406	WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF
407	WRITE(numout,*) ' "ANDREAS" algorithm (Andreas et al. 2015) ln_ANDREAS = ', ln_ANDREAS
408	WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt
409	WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
410	WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
411	WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
412	WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
413	WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk
414	IF(ln_crt_fbk) THEN
415	WRITE(numout,) ' Renault et al. 2020, eq. 10: Stau = Alpha Wnd + Beta'
416	WRITE(numout,*) ' Alpha rn_stau_a = ', rn_stau_a
417	WRITE(numout,*) ' Beta rn_stau_b = ', rn_stau_b
418	ENDIF
419	!
420	WRITE(numout,*)
421	SELECT CASE( nblk ) !* Print the choice of bulk algorithm
422	CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)'
423	CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)'
424	CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)'
425	CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)'
426	CASE( np_ANDREAS ) ; WRITE(numout,*) ' ==>>> "ANDREAS" algorithm (Andreas et al. 2015)'
427	END SELECT
428	!
429	WRITE(numout,*)
430	WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs
431	WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl
432	!
433	WRITE(numout,*)
434	SELECT CASE( nhumi ) !* Print the choice of air humidity
435	CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]'
436	CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]'
437	CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]'
438	END SELECT
439	!
440	!#LB:
441	#if defined key_si3
442	IF( nn_ice > 0 ) THEN
443	WRITE(numout,*)
444	WRITE(numout,*) ' use constant ice-atm bulk transfer coeff. ln_Cx_ice_cst = ', ln_Cx_ice_cst
445	WRITE(numout,*) ' use ice-atm bulk coeff. from Andreas et al., 2005 ln_Cx_ice_AN05 = ', ln_Cx_ice_AN05
446	WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes et al., 2012 ln_Cx_ice_LU12 = ', ln_Cx_ice_LU12
447	WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes & Gryanik, 2015 ln_Cx_ice_LG15 = ', ln_Cx_ice_LG15
448	ENDIF
449	WRITE(numout,*)
450	SELECT CASE( nblk_ice ) !* Print the choice of bulk algorithm
451	CASE( np_ice_cst )
452	WRITE(numout,*) ' ==>>> Constant bulk transfer coefficients over sea-ice:'
453	WRITE(numout,*) ' => Cd_ice, Ce_ice, Ch_ice =', REAL(rn_Cd_i,4), REAL(rn_Ce_i,4), REAL(rn_Ch_i,4)
454	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) ) &
455	& CALL ctl_stop( 'Be realistic in your pick of Cd_ice, Ce_ice & Ch_ice ! (0 < Cx < 1.E-2)')
456	CASE( np_ice_an05 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Andreas et al, 2005'
457	CASE( np_ice_lu12 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes et al, 2012'
458	CASE( np_ice_lg15 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes & Gryanik, 2015'
459	END SELECT
460	#endif
461	!#LB.
462	!
463	ENDIF
464	!
465	END SUBROUTINE sbc_blk_init
466
467
468	SUBROUTINE sbc_blk( kt )
469	!!---------------------------------------------------------------------
470	!! * ROUTINE sbc_blk *
471	!!
472	!! ** Purpose : provide at each time step the surface ocean fluxes
473	!! (momentum, heat, freshwater and runoff)
474	!!
475	!! ** Method :
476	!! (1) READ each fluxes in NetCDF files:
477	!! the wind velocity (i-component) at z=rn_zu (m/s) at T-point
478	!! the wind velocity (j-component) at z=rn_zu (m/s) at T-point
479	!! the specific humidity at z=rn_zqt (kg/kg)
480	!! the air temperature at z=rn_zqt (Kelvin)
481	!! the solar heat (W/m2)
482	!! the Long wave (W/m2)
483	!! the total precipitation (rain+snow) (Kg/m2/s)
484	!! the snow (solid precipitation) (kg/m2/s)
485	!! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds)
486	!! (2) CALL blk_oce_1 and blk_oce_2
487	!!
488	!! C A U T I O N : never mask the surface stress fields
489	!! the stress is assumed to be in the (i,j) mesh referential
490	!!
491	!! ** Action : defined at each time-step at the air-sea interface
492	!! - utau, vtau i- and j-component of the wind stress
493	!! - taum wind stress module at T-point
494	!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
495	!! - qns, qsr non-solar and solar heat fluxes
496	!! - emp upward mass flux (evapo. - precip.)
497	!! - sfx salt flux due to freezing/melting (non-zero only if ice is present)
498	!!
499	!! ** References : Large & Yeager, 2004 / Large & Yeager, 2008
500	!! Brodeau et al. Ocean Modelling 2010
501	!!----------------------------------------------------------------------
502	INTEGER, INTENT(in) :: kt ! ocean time step
503	!!----------------------------------------------------------------------
504	REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp
505	REAL(wp) :: ztmp
506	!!----------------------------------------------------------------------
507	!
508	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
509
510	! Sanity/consistence test on humidity at first time step to detect potential screw-up:
511	IF( kt == nit000 ) THEN
512	IF(lwp) WRITE(numout,*) ''
513	#if defined key_agrif
514	IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ==='
515	#else
516	ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain
517	IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points!
518	ztmp = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztmp ! mean humidity over ocean on proc
519	SELECT CASE( nhumi )
520	CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!)
521	IF( (ztmp < 0._wp) .OR. (ztmp > 0.065) ) ztmp = -1._wp
522	CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K]
523	IF( (ztmp < 110._wp).OR.(ztmp > 320._wp) ) ztmp = -1._wp
524	CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%]
525	IF( (ztmp < 0._wp) .OR.(ztmp > 100._wp) ) ztmp = -1._wp
526	END SELECT
527	IF(ztmp < 0._wp) THEN
528	IF (lwp) WRITE(numout,'(" Mean humidity value found on proc #",i6.6," is: ",f10.5)') narea, ztmp
529	CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', &
530	& ' ==> check the unit in your input files' , &
531	& ' ==> check consistence of namelist choice: specific? relative? dew-point?', &
532	& ' ==> ln_humi_sph -> [kg/kg] \| ln_humi_rlh -> [%] \| ln_humi_dpt -> [K] !!!' )
533	END IF
534	END IF
535	IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ==='
536	#endif
537	IF(lwp) WRITE(numout,*) ''
538	END IF !IF( kt == nit000 )
539	! ! compute the surface ocean fluxes using bulk formulea
540	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
541
542	! Specific humidity of air at z=rn_zqt !
543	SELECT CASE( nhumi )
544	CASE( np_humi_sph )
545	q_air_zt(:,:) = sf(jp_humi )%fnow(:,:,1) ! what we read in file is already a spec. humidity!
546	CASE( np_humi_dpt )
547	IF((kt==nit000).AND.lwp) WRITE(numout,) ' ** sbc_blk() => computing q_air out of dew-point and P !'
548	q_air_zt(:,:) = q_sat( sf(jp_humi )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) )
549	CASE( np_humi_rlh )
550	IF((kt==nit000).AND.lwp) WRITE(numout,) ' ** sbc_blk() => computing q_air out of RH, t_air and slp !' !LBrm
551	q_air_zt(:,:) = q_air_rh( 0.01_wp*sf(jp_humi )%fnow(:,:,1), &
552	& sf(jp_tair )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) !#LB: 0.01 => RH is % percent in file
553	END SELECT
554
555	! POTENTIAL temperature of air at z=rn_zqt
556	! based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
557	! (most reanalysis products provide absolute temp., not potential temp.)
558	IF( ln_tair_pot ) THEN
559	! temperature read into file is already potential temperature, do nothing...
560	theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)
561	ELSE
562	! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...)
563	IF((kt==nit000).AND.lwp) WRITE(numout,) ' ** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!'
564	theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt
565	ENDIF
566	!
567	CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1), & ! <<= in
568	& theta_air_zt(:,:), q_air_zt(:,:), & ! <<= in
569	& sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in
570	& sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in
571	& sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs)
572	& tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out
573
574	CALL blk_oce_2( theta_air_zt(:,:), & ! <<= in
575	& sf(jp_qlw )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1), & ! <<= in
576	& sf(jp_snow )%fnow(:,:,1), tsk_m, & ! <<= in
577	& zsen, zlat, zevp ) ! <=> in out
578	ENDIF
579	!
580	#if defined key_cice
581	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
582	qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1)
583	IF( ln_dm2dc ) THEN
584	qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
585	ELSE
586	qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
587	ENDIF
588	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature instead ????
589	!tatm_ice(:,:) = theta_air_zt(:,:) !#LB: THIS! ?
590
591	qatm_ice(:,:) = q_air_zt(:,:) !#LB:
592
593	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
594	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
595	wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
596	wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
597	ENDIF
598	#endif
599	!
600	END SUBROUTINE sbc_blk
601
602
603	SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! inp
604	& pslp , pst , pu , pv, & ! inp
605	& puatm, pvatm, pdqsr , pdqlw , & ! inp
606	& ptsk , pssq , pcd_du, psen, plat, pevp ) ! out
607	!!---------------------------------------------------------------------
608	!! * ROUTINE blk_oce_1 *
609	!!
610	!! ** Purpose : if ln_blk=T, computes surface momentum, heat and freshwater fluxes
611	!! if ln_abl=T, computes Cd x \|U\|, Ch x \|U\|, Ce x \|U\| for ABL integration
612	!!
613	!! ** Method : bulk formulae using atmospheric fields from :
614	!! if ln_blk=T, atmospheric fields read in sbc_read
615	!! if ln_abl=T, the ABL model at previous time-step
616	!!
617	!! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg)
618	!! - pcd_du : Cd x \|dU\| at T-points (m/s)
619	!! - psen : sensible heat flux (W/m^2)
620	!! - plat : latent heat flux (W/m^2)
621	!! - pevp : evaporation (mm/s) #lolo
622	!!---------------------------------------------------------------------
623	INTEGER , INTENT(in ) :: kt ! time step index
624	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s]
625	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s]
626	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqair ! specific humidity at T-points [kg/kg]
627	REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin]
628	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa]
629	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celsius]
630	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s]
631	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s]
632	REAL(wp), INTENT(in ), DIMENSION(:,:) :: puatm ! surface current seen by the atm at T-point (i-component) [m/s]
633	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pvatm ! surface current seen by the atm at T-point (j-component) [m/s]
634	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqsr ! downwelling solar (shortwave) radiation at surface [W/m^2]
635	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2]
636	REAL(wp), INTENT( out), DIMENSION(:,:) :: ptsk ! skin temp. (or SST if CS & WL not used) [Celsius]
637	REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg]
638	REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du
639	REAL(wp), INTENT( out), DIMENSION(:,:) :: psen
640	REAL(wp), INTENT( out), DIMENSION(:,:) :: plat
641	REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp
642	!
643	INTEGER :: ji, jj ! dummy loop indices
644	REAL(wp) :: zztmp ! local variable
645	REAL(wp) :: zstmax, zstau
646	#if defined key_cyclone
647	REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
648	#endif
649	REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point
650	REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
651	REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean
652	REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean
653	REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean
654	REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2
655	!!---------------------------------------------------------------------
656	!
657	! local scalars ( place there for vector optimisation purposes)
658	! ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K)
659	ptsk(:,:) = pst(:,:) + rt0 ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!)
660
661	! --- cloud cover --- !
662	cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
663
664	! ----------------------------------------------------------------------------- !
665	! 0 Wind components and module at T-point relative to the moving ocean !
666	! ----------------------------------------------------------------------------- !
667
668	! ... components ( U10m - U_oce ) at T-point (unmasked)
669	#if defined key_cyclone
670	zwnd_i(:,:) = 0._wp
671	zwnd_j(:,:) = 0._wp
672	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
673	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
674	zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
675	zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
676	! ... scalar wind at T-point (not masked)
677	wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) )
678	END_2D
679	#else
680	! ... scalar wind module at T-point (not masked)
681	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
682	wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
683	END_2D
684	#endif
685	! ----------------------------------------------------------------------------- !
686	! I Solar FLUX !
687	! ----------------------------------------------------------------------------- !
688
689	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
690	zztmp = 1. - albo
691	IF( ln_dm2dc ) THEN
692	qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
693	ELSE
694	qsr(:,:) = zztmp * pdqsr(:,:) * tmask(:,:,1)
695	ENDIF
696
697
698	! ----------------------------------------------------------------------------- !
699	! II Turbulent FLUXES !
700	! ----------------------------------------------------------------------------- !
701
702	! specific humidity at SST
703	pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
704
705	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
706	!! Backup "bulk SST" and associated spec. hum.
707	zztmp1(:,:) = ptsk(:,:)
708	zztmp2(:,:) = pssq(:,:)
709	ENDIF
710
711	!! Time to call the user-selected bulk parameterization for
712	!! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more...
713	SELECT CASE( nblk )
714
715	CASE( np_NCAR )
716	CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
717	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
718	& nb_iter=nn_iter_algo )
719	!
720	CASE( np_COARE_3p0 )
721	CALL turb_coare3p0( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
722	& ln_skin_cs, ln_skin_wl, &
723	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
724	& nb_iter=nn_iter_algo, &
725	& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
726	!
727	CASE( np_COARE_3p6 )
728	CALL turb_coare3p6( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
729	& ln_skin_cs, ln_skin_wl, &
730	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
731	& nb_iter=nn_iter_algo, &
732	& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
733	!
734	CASE( np_ECMWF )
735	CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
736	& ln_skin_cs, ln_skin_wl, &
737	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
738	& nb_iter=nn_iter_algo, &
739	& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
740	!
741	CASE( np_ANDREAS )
742	CALL turb_andreas ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
743	& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
744	& nb_iter=nn_iter_algo )
745	!
746	CASE DEFAULT
747	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' )
748	!
749	END SELECT
750
751	IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1))
752	IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1))
753	IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1))
754	!! LB: mainly here for debugging purpose:
755	IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ptair-rt0) * tmask(:,:,1)) ! potential temperature at z=zt
756	IF( iom_use('q_zt') ) CALL iom_put("q_zt", pqair * tmask(:,:,1)) ! specific humidity "
757	IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (theta_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu
758	IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity "
759	IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0
760	IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu
761
762	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
763	!! ptsk and pssq have been updated!!!
764	!!
765	!! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq:
766	WHERE ( fr_i(:,:) > 0.001_wp )
767	! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
768	ptsk(:,:) = zztmp1(:,:)
769	pssq(:,:) = zztmp2(:,:)
770	END WHERE
771	END IF
772
773	! Turbulent fluxes over ocean => BULK_FORMULA @ sbc_phy.F90
774	! -------------------------------------------------------------
775
776	IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp
777	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
778	zztmp = zU_zu(ji,jj)
779	wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod
780	pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
781	psen(ji,jj) = zztmp * zch_oce(ji,jj)
782	pevp(ji,jj) = zztmp * zce_oce(ji,jj)
783	rhoa(ji,jj) = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) )
784	END_2D
785	ELSE !== BLK formulation ==! turbulent fluxes computation
786	CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
787	& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
788	& wndm(:,:), zU_zu(:,:), pslp(:,:), &
789	& taum(:,:), psen(:,:), plat(:,:), &
790	& pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac )
791
792	psen(:,:) = psen(:,:) * tmask(:,:,1)
793	plat(:,:) = plat(:,:) * tmask(:,:,1)
794	taum(:,:) = taum(:,:) * tmask(:,:,1)
795	pevp(:,:) = pevp(:,:) * tmask(:,:,1)
796
797	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
798	IF( wndm(ji,jj) > 0._wp ) THEN
799	zztmp = taum(ji,jj) / wndm(ji,jj)
800	#if defined key_cyclone
801	ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
802	ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
803	#else
804	ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
805	ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
806	#endif
807	ELSE
808	ztau_i(ji,jj) = 0._wp
809	ztau_j(ji,jj) = 0._wp
810	ENDIF
811	END_2D
812
813	IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
814	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)
815	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
816	zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) ! stau (<0) must be smaller than zstmax
817	ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) )
818	ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
819	taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) )
820	END_2D
821	ENDIF
822
823	! ... utau, vtau at U- and V_points, resp.
824	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
825	! Note that coastal wind stress is not used in the code... so this extra care has no effect
826	DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T
827	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj ) ) &
828	& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
829	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji ,jj+1) ) &
830	& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
831	END_2D
832
833	IF( ln_crt_fbk ) THEN
834	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. )
835	ELSE
836	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
837	ENDIF
838
839	! Saving open-ocean wind-stress (module and components) on T-points:
840	CALL iom_put( "taum_oce", taum(:,:)*tmask(:,:,1) ) ! output wind stress module
841	!#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])
842	CALL iom_put( "utau_oce", ztau_i(:,:)*tmask(:,:,1) ) ! utau at T-points!
843	CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) ) ! vtau at T-points!
844
845	IF(sn_cfctl%l_prtctl) THEN
846	CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ')
847	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ')
848	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, &
849	& tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
850	CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ')
851	ENDIF
852	!
853	ENDIF !IF( ln_abl )
854
855	ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius
856
857	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
858	CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius
859	CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference...
860	ENDIF
861	!
862	END SUBROUTINE blk_oce_1
863
864
865	SUBROUTINE blk_oce_2( ptair, pdqlw, pprec, psnow, & ! <<= in
866	& ptsk, psen, plat, pevp ) ! <<= in
867	!!---------------------------------------------------------------------
868	!! * ROUTINE blk_oce_2 *
869	!!
870	!! ** Purpose : finalize the momentum, heat and freshwater fluxes computation
871	!! at the ocean surface at each time step knowing Cd, Ch, Ce and
872	!! atmospheric variables (from ABL or external data)
873	!!
874	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
875	!! - vtau : j-component of the stress at V-point (N/m2)
876	!! - taum : Wind stress module at T-point (N/m2)
877	!! - wndm : Wind speed module at T-point (m/s)
878	!! - qsr : Solar heat flux over the ocean (W/m2)
879	!! - qns : Non Solar heat flux over the ocean (W/m2)
880	!! - emp : evaporation minus precipitation (kg/m2/s)
881	!!---------------------------------------------------------------------
882	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair ! potential temperature of air #LB: confirm!
883	REAL(wp), INTENT(in), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2]
884	REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec
885	REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow
886	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptsk ! SKIN surface temperature [Celsius]
887	REAL(wp), INTENT(in), DIMENSION(:,:) :: psen
888	REAL(wp), INTENT(in), DIMENSION(:,:) :: plat
889	REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp
890	!
891	INTEGER :: ji, jj ! dummy loop indices
892	REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable
893	REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! net long wave radiative heat flux
894	!!---------------------------------------------------------------------
895	!
896	! local scalars ( place there for vector optimisation purposes)
897
898	! ----------------------------------------------------------------------------- !
899	! III Net longwave radiative FLUX !
900	! ----------------------------------------------------------------------------- !
901	!! #LB: now moved after Turbulent fluxes because must use the skin temperature rather than bulk SST
902	!! (ptsk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
903	zqlw(:,:) = qlw_net( pdqlw(:,:), ptsk(:,:)+rt0 )
904
905	! ----------------------------------------------------------------------------- !
906	! IV Total FLUXES !
907	! ----------------------------------------------------------------------------- !
908	!
909	emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.)
910	& - pprec(:,:) * rn_pfac ) * tmask(:,:,1)
911	!
912	qns(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) & ! Downward Non Solar
913	& - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip
914	& - pevp(:,:) * ptsk(:,:) * rcp & ! remove evap heat content at SST
915	& + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair
916	& * ( ptair(:,:) - rt0 ) * rcp &
917	& + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
918	& * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi
919	qns(:,:) = qns(:,:) * tmask(:,:,1)
920	!
921	#if defined key_si3
922	qns_oce(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) ! non solar without emp (only needed by SI3)
923	qsr_oce(:,:) = qsr(:,:)
924	#endif
925	!
926	CALL iom_put( "rho_air" , rhoa*tmask(:,:,1) ) ! output air density [kg/m^3]
927	CALL iom_put( "evap_oce" , pevp ) ! evaporation
928	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
929	CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean
930	CALL iom_put( "qla_oce" , plat ) ! output downward latent heat over the ocean
931	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
932	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
933	CALL iom_put( 'snowpre', sprecip ) ! Snow
934	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
935	!
936	IF ( nn_ice == 0 ) THEN
937	CALL iom_put( "qemp_oce" , qns-zqlw-psen-plat ) ! output downward heat content of E-P over the ocean
938	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
939	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
940	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
941	ENDIF
942	!
943	IF(sn_cfctl%l_prtctl) THEN
944	CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ')
945	CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen : ' )
946	CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat : ' )
947	CALL prt_ctl(tab2d_1=qns , clinfo1=' blk_oce_2: qns : ' )
948	CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ')
949	ENDIF
950	!
951	END SUBROUTINE blk_oce_2
952
953
954	#if defined key_si3
955	!!----------------------------------------------------------------------
956	!! 'key_si3' SI3 sea-ice model
957	!!----------------------------------------------------------------------
958	!! blk_ice_1 : provide the air-ice stress
959	!! blk_ice_2 : provide the heat and mass fluxes at air-ice interface
960	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
961	!!----------------------------------------------------------------------
962
963	SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, pqair, pslp , puice, pvice, ptsui, & ! inputs
964	& putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs
965	!!---------------------------------------------------------------------
966	!! * ROUTINE blk_ice_1 *
967	!!
968	!! ** Purpose : provide the surface boundary condition over sea-ice
969	!!
970	!! ** Method : compute momentum using bulk formulation
971	!! formulea, ice variables and read atmospheric fields.
972	!! NB: ice drag coefficient is assumed to be a constant
973	!!---------------------------------------------------------------------
974	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa]
975	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s]
976	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s]
977	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s]
978	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric wind at T-point [m/s]
979	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s]
980	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! "
981	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K]
982	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk
983	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk
984	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl
985	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl
986	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl
987	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl
988	!
989	INTEGER :: ji, jj ! dummy loop indices
990	REAL(wp) :: zootm_su ! sea-ice surface mean temperature
991	REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
992	REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! temporary array
993	!!---------------------------------------------------------------------
994	!
995	! LB: ptsui is in K !!!
996	!
997	! ------------------------------------------------------------ !
998	! Wind module relative to the moving ice ( U10m - U_ice ) !
999	! ------------------------------------------------------------ !
1000	! C-grid ice dynamics : U & V-points (same as ocean)
1001	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1002	wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
1003	END_2D
1004	!
1005	! Make ice-atm. drag dependent on ice concentration
1006
1007
1008	SELECT CASE( nblk_ice )
1009
1010	CASE( np_ice_cst )
1011	! Constant bulk transfer coefficients over sea-ice:
1012	Cd_ice(:,:) = rn_Cd_i
1013	Ch_ice(:,:) = rn_Ch_i
1014	Ce_ice(:,:) = rn_Ce_i
1015	! no height adjustment, keeping zt values:
1016	theta_zu_i(:,:) = ptair(:,:)
1017	q_zu_i(:,:) = pqair(:,:)
1018
1019	CASE( np_ice_an05 ) ! calculate new drag from Lupkes(2015) equations
1020	ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
1021	CALL turb_ice_an05( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, &
1022	& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
1023	!!
1024	CASE( np_ice_lu12 )
1025	ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
1026	CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
1027	& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
1028	!!
1029	CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations
1030	ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
1031	CALL turb_ice_lg15( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
1032	& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
1033	!!
1034	END SELECT
1035
1036	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') ) &
1037	& ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice !
1038
1039	IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp)
1040	IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp)
1041	IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp)
1042
1043
1044	IF( ln_blk ) THEN
1045	! ---------------------------------------------------- !
1046	! Wind stress relative to nonmoving ice ( U10m ) !
1047	! ---------------------------------------------------- !
1048	! supress moving ice in wind stress computation as we don't know how to do it properly...
1049	DO_2D( 0, 1, 0, 1 ) ! at T point
1050	zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
1051	putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
1052	pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
1053	END_2D
1054
1055	!#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!!
1056	IF(iom_use('taum_ice')) CALL iom_put('taum_ice', SQRT( putauiputaui + pvtauipvtaui )*ztmp )
1057	!#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])
1058	IF(iom_use('utau_ice')) CALL iom_put("utau_ice", putaui*ztmp) ! utau at T-points!
1059	IF(iom_use('vtau_ice')) CALL iom_put("vtau_ice", pvtaui*ztmp) ! vtau at T-points!
1060
1061	!
1062	DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean).
1063	!#LB: QUESTION?? so SI3 expects wind stress vector to be provided at U & V points? Not at T-points ?
1064	! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology
1065	zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) )
1066	zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) )
1067	putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) )
1068	pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) )
1069	END_2D
1070	CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
1071	!
1072	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' &
1073	& , tab2d_2=pvtaui , clinfo2=' pvtaui : ' )
1074	ELSE ! ln_abl
1075	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1076	pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj)
1077	pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
1078	pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
1079	END_2D
1080	!#LB:
1081	pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq !
1082	!#LB.
1083	ENDIF !IF( ln_blk )
1084	!
1085	IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
1086	!
1087	END SUBROUTINE blk_ice_1
1088
1089
1090	SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, pqair, pslp, pdqlw, pprec, psnow )
1091	!!---------------------------------------------------------------------
1092	!! * ROUTINE blk_ice_2 *
1093	!!
1094	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
1095	!!
1096	!! ** Method : compute heat and freshwater exchanged
1097	!! between atmosphere and sea-ice using bulk formulation
1098	!! formulea, ice variables and read atmmospheric fields.
1099	!!
1100	!! caution : the net upward water flux has with mm/day unit
1101	!!---------------------------------------------------------------------
1102	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K]
1103	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
1104	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
1105	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
1106	REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair ! potential temperature of air #LB: okay ???
1107	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqair ! specific humidity of air
1108	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp
1109	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pdqlw
1110	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec
1111	REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow
1112	!!
1113	INTEGER :: ji, jj, jl ! dummy loop indices
1114	REAL(wp) :: zst, zst3, zsq ! local variable
1115	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
1116	REAL(wp) :: zztmp, zzblk, zztmp1, z1_rLsub ! - -
1117	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
1118	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
1119	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
1120	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
1121	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
1122	REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2
1123	REAL(wp), DIMENSION(jpi,jpj) :: ztri
1124	!!---------------------------------------------------------------------
1125	!
1126	zcoef_dqlw = 4._wp * emiss_i * stefan ! local scalars
1127	!
1128
1129	zztmp = 1. / ( 1. - albo )
1130	dqla_ice(:,:,:) = 0._wp
1131
1132	! ! ========================== !
1133	DO jl = 1, jpl ! Loop over ice categories !
1134	! ! ========================== !
1135	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1136
1137	zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K]
1138	zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present
1139
1140	! ----------------------------!
1141	! I Radiative FLUXES !
1142	! ----------------------------!
1143	! Short Wave (sw)
1144	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
1145
1146	! Long Wave (lw)
1147	zst3 = zst * zst * zst
1148	z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1)
1149	! lw sensitivity
1150	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
1151
1152	! ----------------------------!
1153	! II Turbulent FLUXES !
1154	! ----------------------------!
1155
1156	! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
1157
1158	! Common term in bulk F. equations...
1159	zzblk = rhoa(ji,jj) * wndm_ice(ji,jj)
1160
1161	! Sensible Heat
1162	zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj)
1163	z_qsb (ji,jj,jl) = zztmp1 * (zst - theta_zu_i(ji,jj))
1164	z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice)
1165
1166	! Latent Heat
1167	zztmp1 = zzblk * rLsub * Ce_ice(ji,jj)
1168	qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ???
1169	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)
1170	! !#LB: dq_sat_dt_ice() in "sbc_phy.F90"
1171	!#LB: without this unjustified "condensation sensure":
1172	!qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj))
1173	!dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
1174
1175
1176
1177	! ----------------------------!
1178	! III Total FLUXES !
1179	! ----------------------------!
1180
1181	! Downward Non Solar flux
1182	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
1183
1184	! Total non solar heat flux sensitivity for ice
1185	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ????
1186
1187	END_2D
1188	!
1189	END DO
1190	!
1191	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]
1192	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]
1193	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
1194	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
1195
1196	! --- evaporation --- !
1197	z1_rLsub = 1._wp / rLsub
1198	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation
1199	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT
1200	zevap (:,:) = emp(:,:) + tprecip(:,:) ! evaporation over ocean !LB: removed rn_efac here, correct???
1201
1202	! --- evaporation minus precipitation --- !
1203	zsnw(:,:) = 0._wp
1204	CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
1205	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
1206	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
1207	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
1208
1209	! --- heat flux associated with emp --- !
1210	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
1211	& + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair
1212	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
1213	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1214	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
1215	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1216
1217	! --- total solar and non solar fluxes --- !
1218	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
1219	& + qemp_ice(:,:) + qemp_oce(:,:)
1220	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
1221
1222	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
1223	qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1224
1225	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
1226	DO jl = 1, jpl
1227	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) rcpi * tmask(:,:,1) )
1228	! ! But we do not have Tice => consider it at 0degC => evap=0
1229	END DO
1230
1231	! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- !
1232	IF( nn_qtrice == 0 ) THEN
1233	! formulation derived from Grenfell and Maykut (1977), where transmission rate
1234	! 1) depends on cloudiness
1235	! 2) is 0 when there is any snow
1236	! 3) tends to 1 for thin ice
1237	ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm
1238	DO jl = 1, jpl
1239	WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
1240	qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
1241	ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm
1242	qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
1243	ELSEWHERE ! zero when hs>0
1244	qtr_ice_top(:,:,jl) = 0._wp
1245	END WHERE
1246	ENDDO
1247	ELSEIF( nn_qtrice == 1 ) THEN
1248	! formulation is derived from the thesis of M. Lebrun (2019).
1249	! It represents the best fit using several sets of observations
1250	! It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
1251	qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:)
1252	ENDIF
1253	!
1254	IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
1255	ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )
1256	IF( iom_use('evap_ao_cea' ) ) CALL iom_put( 'evap_ao_cea' , ztmp(:,:) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
1257	IF( iom_use('hflx_evap_cea') ) CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) ) ! heat flux from evap (cell average)
1258	ENDIF
1259	IF( iom_use('hflx_rain_cea') ) THEN
1260	ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
1261	IF( iom_use('hflx_rain_cea') ) CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) ) ! heat flux from rain (cell average)
1262	ENDIF
1263	IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea') ) THEN
1264	WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 )
1265	ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
1266	ELSEWHERE
1267	ztmp(:,:) = rcp * sst_m(:,:)
1268	ENDWHERE
1269	ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )
1270	IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea' , ztmp2(:,:) ) ! heat flux from snow (cell average)
1271	IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
1272	IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) * zsnw(:,:) ) ! heat flux from snow (over ice)
1273	ENDIF
1274	!
1275	IF(sn_cfctl%l_prtctl) THEN
1276	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
1277	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
1278	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
1279	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
1280	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
1281	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
1282	ENDIF
1283
1284	!#LB:
1285	! air-ice heat flux components that are not written from ice_stp()@icestp.F90:
1286	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
1287	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
1288	IF( iom_use('qlw_ice') ) CALL iom_put( 'qlw_ice', SUM( z_qlw * a_i_b, dim=3 ) )
1289	!#LB.
1290
1291	END SUBROUTINE blk_ice_2
1292
1293
1294	SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
1295	!!---------------------------------------------------------------------
1296	!! * ROUTINE blk_ice_qcn *
1297	!!
1298	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
1299	!! to force sea ice / snow thermodynamics
1300	!! in the case conduction flux is emulated
1301	!!
1302	!! ** Method : compute surface energy balance assuming neglecting heat storage
1303	!! following the 0-layer Semtner (1976) approach
1304	!!
1305	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
1306	!! - qcn_ice : surface inner conduction flux (W/m2)
1307	!!
1308	!!---------------------------------------------------------------------
1309	LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option
1310	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
1311	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
1312	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
1313	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
1314	!
1315	INTEGER , PARAMETER :: nit = 10 ! number of iterations
1316	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
1317	!
1318	INTEGER :: ji, jj, jl ! dummy loop indices
1319	INTEGER :: iter ! local integer
1320	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
1321	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
1322	REAL(wp) :: zqc, zqnet !
1323	REAL(wp) :: zhe, zqa0 !
1324	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
1325	!!---------------------------------------------------------------------
1326
1327	! -------------------------------------!
1328	! I Enhanced conduction factor !
1329	! -------------------------------------!
1330	! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
1331	! Fichefet and Morales Maqueda, JGR 1997
1332	!
1333	zgfac(:,:,:) = 1._wp
1334
1335	IF( ld_virtual_itd ) THEN
1336	!
1337	zfac = 1._wp / ( rn_cnd_s + rcnd_i )
1338	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1339	zfac3 = 2._wp / zepsilon
1340	!
1341	DO jl = 1, jpl
1342	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1343	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness
1344	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
1345	END_2D
1346	END DO
1347	!
1348	ENDIF
1349
1350	! -------------------------------------------------------------!
1351	! II Surface temperature and conduction flux !
1352	! -------------------------------------------------------------!
1353	!
1354	zfac = rcnd_i * rn_cnd_s
1355	!
1356	DO jl = 1, jpl
1357	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1358	!
1359	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
1360	& ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1361	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
1362	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
1363	zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1364	!
1365	DO iter = 1, nit ! --- Iterative loop
1366	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
1367	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
1368	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
1369	END DO
1370	!
1371	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
1372	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1373	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1374	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) ) &
1375	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1376
1377	! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1378	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)
1379
1380	END_2D
1381	!
1382	END DO
1383	!
1384	END SUBROUTINE blk_ice_qcn
1385
1386	#endif
1387
1388	!!======================================================================
1389	END MODULE sbcblk

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

Download in other formats: