Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

sbcblk.F90

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/OCE/SBC/sbcblk.F90

Download in other formats: