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

Last change on this file since 13719 was 13719, checked in by laurent, 4 years ago
Keep up with trunk r13718 + figure generation for STATION_ASF (ocean & ice).
Property svn:keywords set to `Id`
File size: 78.8 KB

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

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

Download in other formats: