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

Last change on this file since 13817 was 13806, checked in by laurent, 4 years ago
Various improvements and cleaning + keep up with "r13801" of the trunk.
Property svn:keywords set to `Id`
File size: 79.1 KB

Rev	Line
[6723]	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	!!=====================================================================
[7163]	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
[6723]	12	!! 3.2 ! 2009-04 (B. Lemaire) Introduce iom_put
	13	!! 3.3 ! 2010-10 (S. Masson) add diurnal cycle
[7163]	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
[12377]	17	!! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
[7163]	18	!! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine
[12377]	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)
[6723]	21	!!----------------------------------------------------------------------
	22
	23	!!----------------------------------------------------------------------
[7163]	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
[12377]	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
[10534]	31	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[6723]	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
[13655]	42	!
[9570]	43	#if defined key_si3
[13655]	44	USE sbc_ice ! Surface boundary condition: ice fields #LB? ok to be in 'key_si3' ???
[13472]	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
[13655]	47	USE sbcblk_algo_ice_lu12
	48	USE sbcblk_algo_ice_lg15
[6723]	49	#endif
[13655]	50	USE sbcblk_algo_ncar ! => turb_ncar : NCAR - (formerly known as CORE, Large & Yeager, 2009)
[12377]	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)
[13655]	54	USE sbcblk_algo_andreas ! => turb_andreas : Andreas et al. 2015
[6727]	55	!
[13655]	56
	57	!
[6723]	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
[13719]	64	USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
[12377]	65
	66
[6723]	67	IMPLICIT NONE
	68	PRIVATE
	69
[7163]	70	PUBLIC sbc_blk_init ! called in sbcmod
	71	PUBLIC sbc_blk ! called in sbcmod
[12377]	72	PUBLIC blk_oce_1 ! called in sbcabl
	73	PUBLIC blk_oce_2 ! called in sbcabl
[9570]	74	#if defined key_si3
[12377]	75	PUBLIC blk_ice_1 ! routine called in icesbc
	76	PUBLIC blk_ice_2 ! routine called in icesbc
[10535]	77	PUBLIC blk_ice_qcn ! routine called in icesbc
[12377]	78	#endif
[6723]	79
[13208]	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
[13472]	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
[6723]	97
[13208]	98	! Warning: keep this structure allocatable for Agrif...
[12377]	99	TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read)
[6723]	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)
[12377]	104	LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013)
	105	LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1)
[13655]	106	LOGICAL :: ln_ANDREAS ! "ANDREAS" algorithm (Andreas et al. 2015)
[6723]	107	!
[13655]	108	!#LB:
[13806]	109	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)
	110	REAL(wp) :: rn_Cd_i, rn_Ce_i, rn_Ch_i ! values for " "
	111	LOGICAL :: ln_Cx_ice_LU12 ! air-ice bulk transfer coefficients based on Lupkes et al., 2012)
	112	LOGICAL :: ln_Cx_ice_LG15 ! air-ice bulk transfer coefficients based on Lupkes & Gryanik, 2015)
[13655]	113	!#LB.
[7355]	114	!
[13208]	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
[13655]	117	REAL(wp) :: rn_stau_b !
[13208]	118	!
[12377]	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	!
[13655]	124	INTEGER :: nn_iter_algo ! Number of iterations in bulk param. algo ("stable ABL + weak wind" requires more)
[6723]	125
[13655]	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
[12377]	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
[13806]	139	LOGICAL :: ln_tair_pot ! temperature read in files ("sn_tair") is already potential temperature (not absolute)
[12377]	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
[6723]	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)
[12377]	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)
[13655]	153	INTEGER, PARAMETER :: np_ANDREAS = 5 ! "ANDREAS" algorithm (Andreas et al. 2015)
[6723]	154
[13655]	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
[6723]	168	!! * Substitutions
[12377]	169	# include "do_loop_substitute.h90"
[6723]	170	!!----------------------------------------------------------------------
[9598]	171	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10069]	172	!! $Id$
[10068]	173	!! Software governed by the CeCILL license (see ./LICENSE)
[6723]	174	!!----------------------------------------------------------------------
	175	CONTAINS
	176
[7355]	177	INTEGER FUNCTION sbc_blk_alloc()
	178	!!-------------------------------------------------------------------
	179	!! * ROUTINE sbc_blk_alloc *
	180	!!-------------------------------------------------------------------
[13655]	181	ALLOCATE( theta_zu(jpi,jpj), q_zu(jpi,jpj), STAT=sbc_blk_alloc )
[10425]	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' )
[7355]	184	END FUNCTION sbc_blk_alloc
[13655]	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
[7355]	197
[9019]	198
[7163]	199	SUBROUTINE sbc_blk_init
	200	!!---------------------------------------------------------------------
	201	!! * ROUTINE sbc_blk_init *
	202	!!
	203	!! ** Purpose : choose and initialize a bulk formulae formulation
	204	!!
[12377]	205	!! ** Method :
[7163]	206	!!
	207	!!----------------------------------------------------------------------
[12377]	208	INTEGER :: jfpr ! dummy loop indice and argument
[7163]	209	INTEGER :: ios, ierror, ioptio ! Local integer
	210	!!
	211	CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files
[13208]	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 ! " "
[13472]	216	TYPE(FLD_N) :: sn_cc, sn_hpgi, sn_hpgj ! " "
[13208]	217	INTEGER :: ipka ! number of levels in the atmospheric variable
[13655]	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, &
[13806]	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_tair_pot, &
	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, &
[13655]	228	& sn_uoatm, sn_voatm, sn_cc, sn_hpgi, sn_hpgj
	229
	230	! cool-skin / warm-layer !LB
[7163]	231	!!---------------------------------------------------------------------
	232	!
[7355]	233	! ! allocate sbc_blk_core array
[13655]	234	IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
[7355]	235	!
[13655]	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	!
[12377]	240	! !** read bulk namelist
[7163]	241	READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
[11536]	242	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' )
[7163]	243	!
	244	READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
[11536]	245	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' )
[7163]	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
[12377]	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
[13655]	264	IF( ln_ANDREAS ) THEN
	265	nblk = np_ANDREAS ; ioptio = ioptio + 1
	266	ENDIF
[7163]	267	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
[12377]	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' )
[13655]	274	IF( ln_ANDREAS ) &
	275	& CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with ANDREAS algorithm' )
[12377]	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' )
[7163]	299	!
	300	IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr
[11536]	301	IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
[12377]	302	IF( sn_qsr%ln_tint ) THEN
[7163]	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
[13655]	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
[7163]	324	! !* set the bulk structure
	325	! !- store namelist information in an array
[12377]	326	!
[13208]	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
[13472]	331	slf_i(jp_slp ) = sn_slp ; slf_i(jp_cc ) = sn_cc
[13208]	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
[7163]	334	!
[13208]	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	!
[7163]	340	! !- allocate the bulk structure
[12377]	341	ALLOCATE( sf(jpfld), STAT=ierror )
[7163]	342	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
[12377]	343	!
[12459]	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	!
[12377]	347	DO jfpr= 1, jpfld
	348	!
[13208]	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)
[13472]	360	IF( jfpr == jp_slp ) THEN
[13208]	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
[13501]	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
[13472]	370	ELSEIF( jfpr == jp_cc ) THEN
	371	sf(jp_cc)%fnow(:,:,1:ipka) = pp_cldf
[13208]	372	ELSE
	373	WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr
	374	CALL ctl_stop( ctmp1 )
	375	ENDIF
[12377]	376	ELSE !-- used field --!
[13208]	377	IF( sf(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) ) ! allocate array for temporal interpolation
[12377]	378	!
[12489]	379	IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 ) &
[13655]	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...' )
[12377]	382	ENDIF
[7163]	383	END DO
	384	!
[12377]	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
[10425]	388	CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
[12377]	389	!drag coefficient read from wave model definable only with mfs bulk formulae and core
	390	ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN
[13655]	391	CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR bulk formulae')
[12377]	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')
[7431]	394	ENDIF
	395	ELSE
[12377]	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
[7431]	403	!
[12377]	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	!
[7163]	412	IF(lwp) THEN !** Control print
	413	!
[12377]	414	WRITE(numout,) ! namelist
[7163]	415	WRITE(numout,*) ' Namelist namsbc_blk (other than data information):'
[13655]	416	WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR
[7163]	417	WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0
[13655]	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
[7163]	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))'
[13208]	426	WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk
	427	IF(ln_crt_fbk) THEN
[13655]	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
[13208]	431	ENDIF
[7163]	432	!
	433	WRITE(numout,*)
	434	SELECT CASE( nblk ) !* Print the choice of bulk algorithm
[9190]	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)'
[12377]	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)'
[13655]	439	CASE( np_ANDREAS ) ; WRITE(numout,*) ' ==>>> "ANDREAS" algorithm (Andreas et al. 2015)'
[7163]	440	END SELECT
	441	!
[12377]	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	!
[13655]	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	!
[7163]	474	ENDIF
	475	!
	476	END SUBROUTINE sbc_blk_init
	477
	478
[6723]	479	SUBROUTINE sbc_blk( kt )
	480	!!---------------------------------------------------------------------
	481	!! * ROUTINE sbc_blk *
	482	!!
	483	!! ** Purpose : provide at each time step the surface ocean fluxes
[9019]	484	!! (momentum, heat, freshwater and runoff)
[6723]	485	!!
[12377]	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
[6723]	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
[12377]	514	!!----------------------------------------------------------------------
[13655]	515	REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp
[12377]	516	REAL(wp) :: ztmp
	517	!!----------------------------------------------------------------------
[6723]	518	!
	519	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
[12377]	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 )
[6723]	550	! ! compute the surface ocean fluxes using bulk formulea
[12377]	551	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
[13655]	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.)
[13806]	569	IF( ln_tair_pot ) THEN
	570	! temperature read into file is already potential temperature, do nothing...
	571	theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)
	572	ELSE
[13655]	573	! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...)
[13806]	574	IF((kt==nit000).AND.lwp) WRITE(numout,) ' ** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!'
[13655]	575	theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt
	576	ENDIF
	577	!
[13208]	578	CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1), & ! <<= in
[13655]	579	& theta_air_zt(:,:), q_air_zt(:,:), & ! <<= in
[13208]	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)
[13655]	583	& tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out
	584
	585	CALL blk_oce_2( theta_air_zt(:,:), & ! <<= in
[13208]	586	& sf(jp_qlw )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1), & ! <<= in
	587	& sf(jp_snow )%fnow(:,:,1), tsk_m, & ! <<= in
[13655]	588	& zsen, zlat, zevp ) ! <=> in out
[12377]	589	ENDIF
	590	!
[6723]	591	#if defined key_cice
	592	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
[7753]	593	qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1)
[12377]	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
[13655]	599	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature instead ????
	600	!tatm_ice(:,:) = theta_air_zt(:,:) !#LB: THIS! ?
[12377]	601
[13655]	602	qatm_ice(:,:) = q_air_zt(:,:) !#LB:
[12377]	603
[7753]	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)
[6723]	608	ENDIF
	609	#endif
	610	!
	611	END SUBROUTINE sbc_blk
	612
	613
[13655]	614	SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! inp
[13208]	615	& pslp , pst , pu , pv, & ! inp
[13655]	616	& puatm, pvatm, pdqsr , pdqlw , & ! inp
	617	& ptsk , pssq , pcd_du, psen, plat, pevp ) ! out
[6723]	618	!!---------------------------------------------------------------------
[12377]	619	!! * ROUTINE blk_oce_1 *
[6723]	620	!!
[12377]	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
[6723]	623	!!
[12377]	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
[6723]	627	!!
[12377]	628	!! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg)
	629	!! - pcd_du : Cd x \|dU\| at T-points (m/s)
[13655]	630	!! - psen : sensible heat flux (W/m^2)
	631	!! - plat : latent heat flux (W/m^2)
	632	!! - pevp : evaporation (mm/s) #lolo
[6723]	633	!!---------------------------------------------------------------------
[12377]	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]
[13655]	637	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqair ! specific humidity at T-points [kg/kg]
[12377]	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]
[13208]	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]
[13655]	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]
[12377]	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]
[13655]	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
[6723]	653	!
	654	INTEGER :: ji, jj ! dummy loop indices
	655	REAL(wp) :: zztmp ! local variable
[13208]	656	REAL(wp) :: zstmax, zstau
	657	#if defined key_cyclone
[9019]	658	REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
[13208]	659	#endif
	660	REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point
[9019]	661	REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
[12377]	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
[6723]	666	!!---------------------------------------------------------------------
	667	!
[7753]	668	! local scalars ( place there for vector optimisation purposes)
[12377]	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!)
[6723]	671
[13472]	672	! --- cloud cover --- !
	673	cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
	674
[6723]	675	! ----------------------------------------------------------------------------- !
	676	! 0 Wind components and module at T-point relative to the moving ocean !
	677	! ----------------------------------------------------------------------------- !
	678
[7753]	679	! ... components ( U10m - U_oce ) at T-point (unmasked)
[12377]	680	#if defined key_cyclone
[7753]	681	zwnd_i(:,:) = 0._wp
	682	zwnd_j(:,:) = 0._wp
[6723]	683	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
[13305]	684	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13806]	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) )
[12377]	689	END_2D
[13208]	690	#else
	691	! ... scalar wind module at T-point (not masked)
[13305]	692	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13806]	693	wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
[13208]	694	END_2D
[6723]	695	#endif
	696	! ----------------------------------------------------------------------------- !
[12377]	697	! I Solar FLUX !
[6723]	698	! ----------------------------------------------------------------------------- !
	699
	700	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
	701	zztmp = 1. - albo
[12377]	702	IF( ln_dm2dc ) THEN
[13655]	703	qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
[12377]	704	ELSE
[13655]	705	qsr(:,:) = zztmp * pdqsr(:,:) * tmask(:,:,1)
[6723]	706	ENDIF
	707
	708
	709	! ----------------------------------------------------------------------------- !
[12377]	710	! II Turbulent FLUXES !
[6723]	711	! ----------------------------------------------------------------------------- !
	712
[12377]	713	! specific humidity at SST
	714	pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
[6723]	715
[12377]	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 )
[13655]	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	!
[12377]	731	CASE( np_COARE_3p0 )
[13655]	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	!
[12377]	738	CASE( np_COARE_3p6 )
[13655]	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	!
[12377]	745	CASE( np_ECMWF )
[13655]	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	!
[6723]	757	CASE DEFAULT
[13655]	758	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' )
	759	!
	760	END SELECT
[12377]	761
[13132]	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:
[13655]	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
[13132]	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
[13655]	772
[12377]	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
[6723]	782	END IF
	783
[13655]	784	! Turbulent fluxes over ocean => BULK_FORMULA @ sbc_phy.F90
[12377]	785	! -------------------------------------------------------------
[6723]	786
[12377]	787	IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp
[13305]	788	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13806]	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) )
[12377]	795	END_2D
	796	ELSE !== BLK formulation ==! turbulent fluxes computation
[13655]	797	CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
[12615]	798	& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
	799	& wndm(:,:), zU_zu(:,:), pslp(:,:), &
[13655]	800	& taum(:,:), psen(:,:), plat(:,:), &
[12615]	801	& pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac )
[12377]	802
	803	psen(:,:) = psen(:,:) * tmask(:,:,1)
[13655]	804	plat(:,:) = plat(:,:) * tmask(:,:,1)
[12377]	805	taum(:,:) = taum(:,:) * tmask(:,:,1)
	806	pevp(:,:) = pevp(:,:) * tmask(:,:,1)
	807
[13305]	808	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13655]	809	IF( wndm(ji,jj) > 0._wp ) THEN
	810	zztmp = taum(ji,jj) / wndm(ji,jj)
[13208]	811	#if defined key_cyclone
[13655]	812	ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
	813	ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
[13208]	814	#else
[13655]	815	ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
	816	ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
[13208]	817	#endif
[13655]	818	ELSE
	819	ztau_i(ji,jj) = 0._wp
	820	ztau_j(ji,jj) = 0._wp
	821	ENDIF
[13208]	822	END_2D
[12377]	823
[13208]	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)
[13295]	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
[13806]	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) )
[13208]	831	END_2D
	832	ENDIF
[12377]	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
[12925]	836	! Note that coastal wind stress is not used in the code... so this extra care has no effect
[13295]	837	DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T
[13806]	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))
[12377]	842	END_2D
[6723]	843
[13208]	844	IF( ln_crt_fbk ) THEN
	845	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. )
	846	ELSE
	847	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
	848	ENDIF
	849
[13655]	850	CALL iom_put( "taum_oce", taum*tmask(:,:,1) ) ! output wind stress module
[13208]	851
[12377]	852	IF(sn_cfctl%l_prtctl) THEN
[13655]	853	CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ')
	854	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ')
	855	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, &
	856	& tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
	857	CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ')
[12377]	858	ENDIF
	859	!
	860	ENDIF !IF( ln_abl )
[13655]	861
[12377]	862	ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius
[13655]	863
[12377]	864	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
	865	CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius
	866	CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference...
	867	ENDIF
	868	!
	869	END SUBROUTINE blk_oce_1
[6723]	870
[13806]	871
[13655]	872	SUBROUTINE blk_oce_2( ptair, pdqlw, pprec, psnow, & ! <<= in
	873	& ptsk, psen, plat, pevp ) ! <<= in
[12377]	874	!!---------------------------------------------------------------------
	875	!! * ROUTINE blk_oce_2 *
	876	!!
	877	!! ** Purpose : finalize the momentum, heat and freshwater fluxes computation
	878	!! at the ocean surface at each time step knowing Cd, Ch, Ce and
	879	!! atmospheric variables (from ABL or external data)
	880	!!
	881	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
	882	!! - vtau : j-component of the stress at V-point (N/m2)
	883	!! - taum : Wind stress module at T-point (N/m2)
	884	!! - wndm : Wind speed module at T-point (m/s)
	885	!! - qsr : Solar heat flux over the ocean (W/m2)
	886	!! - qns : Non Solar heat flux over the ocean (W/m2)
	887	!! - emp : evaporation minus precipitation (kg/m2/s)
	888	!!---------------------------------------------------------------------
[13655]	889	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair ! potential temperature of air #LB: confirm!
	890	REAL(wp), INTENT(in), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2]
[12377]	891	REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec
	892	REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow
	893	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptsk ! SKIN surface temperature [Celsius]
	894	REAL(wp), INTENT(in), DIMENSION(:,:) :: psen
[13655]	895	REAL(wp), INTENT(in), DIMENSION(:,:) :: plat
[12377]	896	REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp
	897	!
	898	INTEGER :: ji, jj ! dummy loop indices
	899	REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable
[13655]	900	REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! net long wave radiative heat flux
[12377]	901	!!---------------------------------------------------------------------
	902	!
	903	! local scalars ( place there for vector optimisation purposes)
[6723]	904
[12377]	905	! ----------------------------------------------------------------------------- !
	906	! III Net longwave radiative FLUX !
	907	! ----------------------------------------------------------------------------- !
[13655]	908	!! #LB: now moved after Turbulent fluxes because must use the skin temperature rather than bulk SST
	909	!! (ptsk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
	910	zqlw(:,:) = qlw_net( pdqlw(:,:), ptsk(:,:)+rt0 )
[12377]	911
[6723]	912	! ----------------------------------------------------------------------------- !
[12377]	913	! IV Total FLUXES !
[6723]	914	! ----------------------------------------------------------------------------- !
	915	!
[12377]	916	emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.)
	917	& - pprec(:,:) * rn_pfac ) * tmask(:,:,1)
[7753]	918	!
[13655]	919	qns(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) & ! Downward Non Solar
[12377]	920	& - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip
	921	& - pevp(:,:) * ptsk(:,:) * rcp & ! remove evap heat content at SST
	922	& + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair
	923	& * ( ptair(:,:) - rt0 ) * rcp &
	924	& + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
	925	& * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi
[9727]	926	qns(:,:) = qns(:,:) * tmask(:,:,1)
[7753]	927	!
[9570]	928	#if defined key_si3
[13655]	929	qns_oce(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) ! non solar without emp (only needed by SI3)
[7753]	930	qsr_oce(:,:) = qsr(:,:)
[6723]	931	#endif
	932	!
[12377]	933	CALL iom_put( "rho_air" , rhoa*tmask(:,:,1) ) ! output air density [kg/m^3]
	934	CALL iom_put( "evap_oce" , pevp ) ! evaporation
	935	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
	936	CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean
[13655]	937	CALL iom_put( "qla_oce" , plat ) ! output downward latent heat over the ocean
[12377]	938	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
	939	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
	940	CALL iom_put( 'snowpre', sprecip ) ! Snow
	941	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
	942	!
[6723]	943	IF ( nn_ice == 0 ) THEN
[13655]	944	CALL iom_put( "qemp_oce" , qns-zqlw-psen-plat ) ! output downward heat content of E-P over the ocean
[12377]	945	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
	946	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
	947	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
[6723]	948	ENDIF
	949	!
[12377]	950	IF(sn_cfctl%l_prtctl) THEN
	951	CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ')
[13655]	952	CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen : ' )
	953	CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat : ' )
	954	CALL prt_ctl(tab2d_1=qns , clinfo1=' blk_oce_2: qns : ' )
[12377]	955	CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ')
[6723]	956	ENDIF
	957	!
[12377]	958	END SUBROUTINE blk_oce_2
[6723]	959
	960
[9570]	961	#if defined key_si3
[9019]	962	!!----------------------------------------------------------------------
[9570]	963	!! 'key_si3' SI3 sea-ice model
[9019]	964	!!----------------------------------------------------------------------
[12377]	965	!! blk_ice_1 : provide the air-ice stress
	966	!! blk_ice_2 : provide the heat and mass fluxes at air-ice interface
[10534]	967	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[9019]	968	!!----------------------------------------------------------------------
	969
[13655]	970	SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, pqair, pslp , puice, pvice, ptsui, & ! inputs
[12377]	971	& putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs
[6723]	972	!!---------------------------------------------------------------------
[12377]	973	!! * ROUTINE blk_ice_1 *
[6723]	974	!!
	975	!! ** Purpose : provide the surface boundary condition over sea-ice
	976	!!
	977	!! ** Method : compute momentum using bulk formulation
	978	!! formulea, ice variables and read atmospheric fields.
	979	!! NB: ice drag coefficient is assumed to be a constant
	980	!!---------------------------------------------------------------------
[12377]	981	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa]
	982	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s]
	983	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s]
	984	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s]
[13655]	985	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric wind at T-point [m/s]
[12377]	986	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s]
	987	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! "
	988	REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K]
	989	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk
	990	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk
	991	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl
	992	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl
	993	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl
	994	REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl
	995	!
[6723]	996	INTEGER :: ji, jj ! dummy loop indices
[12377]	997	REAL(wp) :: zootm_su ! sea-ice surface mean temperature
[13655]	998	REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
	999	REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! temporary array
[6723]	1000	!!---------------------------------------------------------------------
	1001	!
[13655]	1002	! LB: ptsui is in K !!!
	1003	!
[9019]	1004	! ------------------------------------------------------------ !
	1005	! Wind module relative to the moving ice ( U10m - U_ice ) !
	1006	! ------------------------------------------------------------ !
[9767]	1007	! C-grid ice dynamics : U & V-points (same as ocean)
[13305]	1008	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13655]	1009	wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
[12377]	1010	END_2D
[9767]	1011	!
[9019]	1012	! Make ice-atm. drag dependent on ice concentration
[6723]	1013
[13655]	1014
	1015	SELECT CASE( nblk_ice )
	1016
	1017	CASE( np_ice_cst )
	1018	! Constant bulk transfer coefficients over sea-ice:
	1019	Cd_ice(:,:) = rn_Cd_i
	1020	Ch_ice(:,:) = rn_Ch_i
	1021	Ce_ice(:,:) = rn_Ce_i
	1022	! no height adjustment, keeping zt values:
	1023	theta_zu_i(:,:) = ptair(:,:)
	1024	q_zu_i(:,:) = pqair(:,:)
	1025
	1026	CASE( np_ice_lu12 )
	1027	ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
	1028	CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
	1029	& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
	1030	!!
	1031	CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations
	1032	ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
	1033	CALL turb_ice_lg15( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
	1034	& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
	1035	!!
	1036	END SELECT
	1037
	1038	IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ai') ) &
	1039	& ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice !
	1040
	1041	IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp)
	1042	IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp)
	1043	IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp)
	1044
	1045
[12377]	1046	IF( ln_blk ) THEN
[13208]	1047	! ---------------------------------------------------- !
	1048	! Wind stress relative to nonmoving ice ( U10m ) !
	1049	! ---------------------------------------------------- !
	1050	! supress moving ice in wind stress computation as we don't know how to do it properly...
[13655]	1051	DO_2D( 0, 1, 0, 1 ) ! at T point
	1052	zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
	1053	putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
	1054	pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
[12377]	1055	END_2D
[13655]	1056	!#LB: saving the module of the ai wind-stress: NOT weighted by the ice concentration !!!
	1057	IF(iom_use('taum_ai')) CALL iom_put( 'taum_ai', SQRT( putauiputaui + pvtauipvtaui )*ztmp )
[12925]	1058	!
[13295]	1059	DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean).
[13655]	1060	!#LB: QUESTION?? so SI3 expects wind stress vector to be provided at U & V points? Not at T-points ?
	1061	! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology
	1062	zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) )
	1063	zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) )
	1064	putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) )
	1065	pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) )
[12925]	1066	END_2D
[13226]	1067	CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
[12377]	1068	!
	1069	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' &
	1070	& , tab2d_2=pvtaui , clinfo2=' pvtaui : ' )
[13214]	1071	ELSE ! ln_abl
[13305]	1072	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13655]	1073	pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj)
	1074	pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
	1075	pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
[12377]	1076	END_2D
[13655]	1077	!#LB:
	1078	pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq !
	1079	!#LB.
	1080	ENDIF !IF( ln_blk )
[9019]	1081	!
[12377]	1082	IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
[9767]	1083	!
[12377]	1084	END SUBROUTINE blk_ice_1
[6723]	1085
	1086
[13655]	1087	SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, pqair, pslp, pdqlw, pprec, psnow )
[6723]	1088	!!---------------------------------------------------------------------
[12377]	1089	!! * ROUTINE blk_ice_2 *
[6723]	1090	!!
[9019]	1091	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
[6723]	1092	!!
	1093	!! ** Method : compute heat and freshwater exchanged
	1094	!! between atmosphere and sea-ice using bulk formulation
	1095	!! formulea, ice variables and read atmmospheric fields.
	1096	!!
	1097	!! caution : the net upward water flux has with mm/day unit
	1098	!!---------------------------------------------------------------------
[12377]	1099	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K]
[9019]	1100	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
	1101	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
[6727]	1102	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
[13655]	1103	REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair ! potential temperature of air #LB: okay ???
	1104	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqair ! specific humidity of air
[12377]	1105	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp
[13655]	1106	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pdqlw
[12377]	1107	REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec
	1108	REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow
[6723]	1109	!!
[6727]	1110	INTEGER :: ji, jj, jl ! dummy loop indices
[13655]	1111	REAL(wp) :: zst, zst3, zsq ! local variable
[6727]	1112	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
[13655]	1113	REAL(wp) :: zztmp, zzblk, zztmp1, z1_rLsub ! - -
[9019]	1114	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
	1115	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
	1116	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
	1117	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
[9656]	1118	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
[12276]	1119	REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2
[13472]	1120	REAL(wp), DIMENSION(jpi,jpj) :: ztri
[6723]	1121	!!---------------------------------------------------------------------
	1122	!
[13655]	1123	zcoef_dqlw = 4._wp * emiss_i * stefan ! local scalars
[6723]	1124	!
[13655]	1125
[6723]	1126	zztmp = 1. / ( 1. - albo )
[13655]	1127	dqla_ice(:,:,:) = 0._wp
	1128
[7753]	1129	! ! ========================== !
	1130	DO jl = 1, jpl ! Loop over ice categories !
	1131	! ! ========================== !
[13655]	1132	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[6723]	1133
[13655]	1134	zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K]
	1135	zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present
[6723]	1136
[13655]	1137	! ----------------------------!
	1138	! I Radiative FLUXES !
	1139	! ----------------------------!
	1140	! Short Wave (sw)
	1141	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
[6723]	1142
[13655]	1143	! Long Wave (lw)
	1144	zst3 = zst * zst * zst
	1145	z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1)
	1146	! lw sensitivity
	1147	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
[6723]	1148
[13655]	1149	! ----------------------------!
	1150	! II Turbulent FLUXES !
	1151	! ----------------------------!
	1152
	1153	! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
	1154
	1155	! Common term in bulk F. equations...
	1156	zzblk = rhoa(ji,jj) * wndm_ice(ji,jj)
	1157
	1158	! Sensible Heat
	1159	zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj)
	1160	z_qsb (ji,jj,jl) = zztmp1 * (zst - theta_zu_i(ji,jj))
	1161	z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice)
	1162
	1163	! Latent Heat
	1164	zztmp1 = zzblk * rLsub * Ce_ice(ji,jj)
	1165	qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ???
	1166	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)
	1167	! !#LB: dq_sat_dt_ice() in "sbc_phy.F90"
	1168	!#LB: without this unjustified "condensation sensure":
	1169	!qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj))
	1170	!dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
	1171
	1172
	1173
	1174	! ----------------------------!
	1175	! III Total FLUXES !
	1176	! ----------------------------!
	1177
	1178	! Downward Non Solar flux
	1179	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
	1180
	1181	! Total non solar heat flux sensitivity for ice
	1182	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ????
	1183
	1184	END_2D
[6723]	1185	!
	1186	END DO
	1187	!
[12377]	1188	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]
	1189	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]
	1190	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
	1191	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
[6723]	1192
	1193	! --- evaporation --- !
[9935]	1194	z1_rLsub = 1._wp / rLsub
	1195	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation
	1196	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT
[12629]	1197	zevap (:,:) = emp(:,:) + tprecip(:,:) ! evaporation over ocean !LB: removed rn_efac here, correct???
[6723]	1198
[7753]	1199	! --- evaporation minus precipitation --- !
	1200	zsnw(:,:) = 0._wp
[13472]	1201	CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
[9019]	1202	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
[7753]	1203	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
	1204	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
[6723]	1205
[7753]	1206	! --- heat flux associated with emp --- !
[9019]	1207	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
[12377]	1208	& + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair
[7753]	1209	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
[12377]	1210	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[7753]	1211	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
[12377]	1212	& ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	1213
[7753]	1214	! --- total solar and non solar fluxes --- !
[9019]	1215	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
	1216	& + qemp_ice(:,:) + qemp_oce(:,:)
	1217	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
[6723]	1218
[7753]	1219	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
[12377]	1220	qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	1221
[7504]	1222	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
	1223	DO jl = 1, jpl
[9935]	1224	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) rcpi * tmask(:,:,1) )
[12377]	1225	! ! But we do not have Tice => consider it at 0degC => evap=0
[7504]	1226	END DO
	1227
[13472]	1228	! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- !
	1229	IF( nn_qtrice == 0 ) THEN
	1230	! formulation derived from Grenfell and Maykut (1977), where transmission rate
	1231	! 1) depends on cloudiness
	1232	! 2) is 0 when there is any snow
	1233	! 3) tends to 1 for thin ice
	1234	ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm
	1235	DO jl = 1, jpl
[13655]	1236	WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
[13472]	1237	qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
	1238	ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm
	1239	qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
	1240	ELSEWHERE ! zero when hs>0
[13655]	1241	qtr_ice_top(:,:,jl) = 0._wp
[13472]	1242	END WHERE
	1243	ENDDO
	1244	ELSEIF( nn_qtrice == 1 ) THEN
	1245	! formulation is derived from the thesis of M. Lebrun (2019).
	1246	! It represents the best fit using several sets of observations
	1247	! It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
	1248	qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:)
	1249	ENDIF
[6723]	1250	!
[12276]	1251	IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
[12377]	1252	ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )
	1253	IF( iom_use('evap_ao_cea' ) ) CALL iom_put( 'evap_ao_cea' , ztmp(:,:) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
	1254	IF( iom_use('hflx_evap_cea') ) CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) ) ! heat flux from evap (cell average)
[12276]	1255	ENDIF
	1256	IF( iom_use('hflx_rain_cea') ) THEN
	1257	ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
[12377]	1258	IF( iom_use('hflx_rain_cea') ) CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) ) ! heat flux from rain (cell average)
[12276]	1259	ENDIF
	1260	IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea') ) THEN
[12377]	1261	WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 )
	1262	ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
	1263	ELSEWHERE
	1264	ztmp(:,:) = rcp * sst_m(:,:)
	1265	ENDWHERE
	1266	ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )
	1267	IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea' , ztmp2(:,:) ) ! heat flux from snow (cell average)
	1268	IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
	1269	IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) * zsnw(:,:) ) ! heat flux from snow (over ice)
[12276]	1270	ENDIF
	1271	!
[12377]	1272	IF(sn_cfctl%l_prtctl) THEN
[6723]	1273	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
	1274	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
	1275	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
	1276	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
	1277	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
	1278	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
	1279	ENDIF
[13655]	1280
	1281	!#LB:
	1282	! air-ice heat flux components that are not written from ice_stp()@icestp.F90:
	1283	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
	1284	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
	1285	IF( iom_use('qlw_ice') ) CALL iom_put( 'qlw_ice', SUM( z_qlw * a_i_b, dim=3 ) )
	1286	!#LB.
	1287
[12377]	1288	END SUBROUTINE blk_ice_2
[6723]	1289
[12377]	1290
[10531]	1291	SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
[9019]	1292	!!---------------------------------------------------------------------
	1293	!! * ROUTINE blk_ice_qcn *
[6723]	1294	!!
[9019]	1295	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
	1296	!! to force sea ice / snow thermodynamics
[10534]	1297	!! in the case conduction flux is emulated
[12377]	1298	!!
[9019]	1299	!! ** Method : compute surface energy balance assuming neglecting heat storage
	1300	!! following the 0-layer Semtner (1976) approach
[6727]	1301	!!
[9019]	1302	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
	1303	!! - qcn_ice : surface inner conduction flux (W/m2)
	1304	!!
	1305	!!---------------------------------------------------------------------
[10531]	1306	LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option
[9076]	1307	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
	1308	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
	1309	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
	1310	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
[6723]	1311	!
[9019]	1312	INTEGER , PARAMETER :: nit = 10 ! number of iterations
	1313	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
[6723]	1314	!
[9019]	1315	INTEGER :: ji, jj, jl ! dummy loop indices
	1316	INTEGER :: iter ! local integer
	1317	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
	1318	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
	1319	REAL(wp) :: zqc, zqnet !
	1320	REAL(wp) :: zhe, zqa0 !
	1321	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
	1322	!!---------------------------------------------------------------------
[12377]	1323
[9019]	1324	! -------------------------------------!
	1325	! I Enhanced conduction factor !
	1326	! -------------------------------------!
[10531]	1327	! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
[9019]	1328	! Fichefet and Morales Maqueda, JGR 1997
[6723]	1329	!
[9019]	1330	zgfac(:,:,:) = 1._wp
[12377]	1331
[10531]	1332	IF( ld_virtual_itd ) THEN
[9019]	1333	!
[9935]	1334	zfac = 1._wp / ( rn_cnd_s + rcnd_i )
[9019]	1335	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
	1336	zfac3 = 2._wp / zepsilon
[12377]	1337	!
	1338	DO jl = 1, jpl
[13305]	1339	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13806]	1340	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness
	1341	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
[12377]	1342	END_2D
[9019]	1343	END DO
[12377]	1344	!
[10531]	1345	ENDIF
[12377]	1346
[9019]	1347	! -------------------------------------------------------------!
	1348	! II Surface temperature and conduction flux !
	1349	! -------------------------------------------------------------!
[6723]	1350	!
[9935]	1351	zfac = rcnd_i * rn_cnd_s
[6723]	1352	!
[9019]	1353	DO jl = 1, jpl
[13305]	1354	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13806]	1355	!
	1356	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
	1357	& ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
	1358	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
	1359	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
	1360	zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
	1361	!
	1362	DO iter = 1, nit ! --- Iterative loop
	1363	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
	1364	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
	1365	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
	1366	END DO
	1367	!
	1368	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
	1369	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
	1370	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
	1371	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) ) &
	1372	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
[6723]	1373
[13806]	1374	! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
	1375	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)
[9938]	1376
[12377]	1377	END_2D
[9019]	1378	!
[12377]	1379	END DO
	1380	!
[9019]	1381	END SUBROUTINE blk_ice_qcn
[6723]	1382
[7355]	1383	#endif
	1384
[6723]	1385	!!======================================================================
	1386	END MODULE sbcblk

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

Download in other formats: