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/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90 @ 11638

Last change on this file since 11638 was 11638, checked in by laurent, 5 years ago
LB: preliminary inclusion of "STATION_ASF" test-case!
Property svn:keywords set to `Id`
File size: 71.5 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
[11111]	17	!! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
[7163]	18	!! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine
[11217]	19	!! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle)
[6723]	20	!!----------------------------------------------------------------------
	21
	22	!!----------------------------------------------------------------------
[7163]	23	!! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition
	24	!! sbc_blk : bulk formulation as ocean surface boundary condition
[6727]	25	!! blk_oce : computes momentum, heat and freshwater fluxes over ocean
[11217]	26	!! sea-ice case only :
[9019]	27	!! blk_ice_tau : provide the air-ice stress
	28	!! blk_ice_flx : provide the heat and mass fluxes at air-ice interface
[10534]	29	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[9019]	30	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
[11217]	31	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
[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 sbc_ice ! Surface boundary condition: ice fields
	42	USE lib_fortran ! to use key_nosignedzero
[9570]	43	#if defined key_si3
[10511]	44	USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
[9019]	45	USE icethd_dh ! for CALL ice_thd_snwblow
[6723]	46	#endif
[11217]	47	USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009)
[11284]	48	USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003)
	49	USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013)
[11615]	50	USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 45r1)
[6727]	51	!
[6723]	52	USE iom ! I/O manager library
	53	USE in_out_manager ! I/O manager
	54	USE lib_mpp ! distribued memory computing library
	55	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	56	USE prtctl ! Print control
	57
[11182]	58	USE sbcblk_phy !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc...
[11111]	59
[11217]	60
[6723]	61	IMPLICIT NONE
	62	PRIVATE
	63
[7163]	64	PUBLIC sbc_blk_init ! called in sbcmod
	65	PUBLIC sbc_blk ! called in sbcmod
[9570]	66	#if defined key_si3
[10535]	67	PUBLIC blk_ice_tau ! routine called in icesbc
	68	PUBLIC blk_ice_flx ! routine called in icesbc
	69	PUBLIC blk_ice_qcn ! routine called in icesbc
[11217]	70	#endif
[6723]	71
	72	INTEGER , PARAMETER :: jpfld =10 ! maximum number of files to read
	73	INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
	74	INTEGER , PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
[7163]	75	INTEGER , PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin)
	76	INTEGER , PARAMETER :: jp_humi = 4 ! index of specific humidity ( % )
	77	INTEGER , PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2)
	78	INTEGER , PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2)
[6723]	79	INTEGER , PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
	80	INTEGER , PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s)
[7421]	81	INTEGER , PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa)
[6723]	82	INTEGER , PARAMETER :: jp_tdif =10 ! index of tau diff associated to HF tau (N/m2) at T-point
	83
	84	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
	85
	86	! !!* Namelist namsbc_blk : bulk parameters
	87	LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008)
	88	LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003)
[11284]	89	LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013)
[11615]	90	LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1)
[6723]	91	!
	92	LOGICAL :: ln_taudif ! logical flag to use the "mean of stress module - module of mean stress" data
	93	REAL(wp) :: rn_pfac ! multiplication factor for precipitation
[9019]	94	REAL(wp) :: rn_efac ! multiplication factor for evaporation
	95	REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress
[6723]	96	REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
	97	REAL(wp) :: rn_zu ! z(u) : height of wind measurements
[11217]	98	!!gm ref namelist initialize it so remove the setting to false below
[9019]	99	LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
	100	LOGICAL :: ln_Cd_L15 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
[7355]	101	!
[9019]	102	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_atm ! transfer coefficient for momentum (tau)
	103	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ch_atm ! transfer coefficient for sensible heat (Q_sens)
	104	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ce_atm ! tansfert coefficient for evaporation (Q_lat)
	105	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_zu ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
	106	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_zu ! air spec. hum. at wind speed height (needed by Lupkes 2015 bulk scheme)
	107	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
[6723]	108
[11111]	109	!LB:
[11615]	110	LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB
	111	LOGICAL :: ln_skin_wl ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
[11284]	112	LOGICAL :: ln_humi_sph ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB
	113	LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB
	114	LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB
	115	!
	116	INTEGER :: nhumi ! choice of the bulk algorithm
	117	! ! associated indices:
	118	INTEGER, PARAMETER :: np_humi_sph = 1
	119	INTEGER, PARAMETER :: np_humi_dpt = 2
	120	INTEGER, PARAMETER :: np_humi_rlh = 3
[11111]	121	!LB.
[11217]	122
[6723]	123	INTEGER :: nblk ! choice of the bulk algorithm
	124	! ! associated indices:
	125	INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008)
	126	INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003)
[11284]	127	INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013)
[11615]	128	INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 45r1)
[6723]	129
	130	!! * Substitutions
	131	# include "vectopt_loop_substitute.h90"
	132	!!----------------------------------------------------------------------
[9598]	133	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10069]	134	!! $Id$
[10068]	135	!! Software governed by the CeCILL license (see ./LICENSE)
[6723]	136	!!----------------------------------------------------------------------
	137	CONTAINS
	138
[7355]	139	INTEGER FUNCTION sbc_blk_alloc()
	140	!!-------------------------------------------------------------------
	141	!! * ROUTINE sbc_blk_alloc *
	142	!!-------------------------------------------------------------------
[9019]	143	ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
	144	& cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
[7355]	145	!
[10425]	146	CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
	147	IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
[7355]	148	END FUNCTION sbc_blk_alloc
	149
[11111]	150	!LB:
	151	INTEGER FUNCTION sbc_blk_cswl_alloc()
	152	!!-------------------------------------------------------------------
	153	!! * ROUTINE sbc_blk_cswl_alloc *
	154	!!-------------------------------------------------------------------
[11637]	155	!WRITE(numout,) '** LB: allocating tsk!'
[11111]	156	ALLOCATE( tsk(jpi,jpj), STAT=sbc_blk_cswl_alloc )
[11637]	157	!WRITE(numout,) '** LB: done!'
[11111]	158	CALL mpp_sum ( 'sbcblk', sbc_blk_cswl_alloc )
	159	IF( sbc_blk_cswl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_cswl_alloc: failed to allocate arrays' )
	160	END FUNCTION sbc_blk_cswl_alloc
	161	!LB.
[9019]	162
[11111]	163
[7163]	164	SUBROUTINE sbc_blk_init
	165	!!---------------------------------------------------------------------
	166	!! * ROUTINE sbc_blk_init *
	167	!!
	168	!! ** Purpose : choose and initialize a bulk formulae formulation
	169	!!
[11217]	170	!! ** Method :
[7163]	171	!!
	172	!!----------------------------------------------------------------------
	173	INTEGER :: ifpr, jfld ! dummy loop indice and argument
	174	INTEGER :: ios, ierror, ioptio ! Local integer
	175	!!
	176	CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files
	177	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
	178	TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read
	179	TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " "
	180	TYPE(FLD_N) :: sn_slp , sn_tdif ! " "
	181	NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields
	182	& sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, &
[11284]	183	& ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm
[11217]	184	& cn_dir , ln_taudif, rn_zqt, rn_zu, &
[11111]	185	& rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, &
[11615]	186	& ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB
[7163]	187	!!---------------------------------------------------------------------
	188	!
[7355]	189	! ! allocate sbc_blk_core array
	190	IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
	191	!
[11217]	192	! !** read bulk namelist
[7163]	193	REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters
	194	READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
[9168]	195	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp )
[7163]	196	!
	197	REWIND( numnam_cfg ) !* Namelist namsbc_blk in configuration namelist : bulk parameters
	198	READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
[9168]	199	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist', lwp )
[7163]	200	!
	201	IF(lwm) WRITE( numond, namsbc_blk )
	202	!
	203	! !** initialization of the chosen bulk formulae (+ check)
	204	! !* select the bulk chosen in the namelist and check the choice
[11217]	205	ioptio = 0
	206	IF( ln_NCAR ) THEN
	207	nblk = np_NCAR ; ioptio = ioptio + 1
	208	ENDIF
	209	IF( ln_COARE_3p0 ) THEN
	210	nblk = np_COARE_3p0 ; ioptio = ioptio + 1
	211	ENDIF
[11284]	212	IF( ln_COARE_3p6 ) THEN
	213	nblk = np_COARE_3p6 ; ioptio = ioptio + 1
[11217]	214	ENDIF
	215	IF( ln_ECMWF ) THEN
	216	nblk = np_ECMWF ; ioptio = ioptio + 1
	217	ENDIF
[11284]	218	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
	219
	220
	221
[11111]	222	!LB:
	223	! !** initialization of the cool-skin / warm-layer parametrization
[11615]	224	IF( ln_skin_cs .OR. ln_skin_wl ) THEN
[11111]	225	IF ( ln_NCAR ) CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' )
	226	! ! allocate array(s) for cool-skin/warm-layer param.
	227	IF( sbc_blk_cswl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
	228	END IF
[11284]	229	!
	230	ioptio = 0
	231	IF( ln_humi_sph ) THEN
	232	nhumi = np_humi_sph ; ioptio = ioptio + 1
	233	ENDIF
	234	IF( ln_humi_dpt ) THEN
	235	nhumi = np_humi_dpt ; ioptio = ioptio + 1
	236	ENDIF
	237	IF( ln_humi_rlh ) THEN
	238	nhumi = np_humi_rlh ; ioptio = ioptio + 1
	239	ENDIF
	240	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' )
[11111]	241	!LB.
[11217]	242
[7163]	243	!
	244	IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr
	245	IF( sn_qsr%nfreqh /= 24 ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
[11217]	246	IF( sn_qsr%ln_tint ) THEN
[7163]	247	CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', &
	248	& ' ==> We force time interpolation = .false. for qsr' )
[11615]	249	sn_qsr%ln_tint = .false.
[7163]	250	ENDIF
	251	ENDIF
	252	! !* set the bulk structure
	253	! !- store namelist information in an array
	254	slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj
	255	slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw
	256	slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi
	257	slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow
	258	slf_i(jp_slp) = sn_slp ; slf_i(jp_tdif) = sn_tdif
	259	!
	260	lhftau = ln_taudif !- add an extra field if HF stress is used
	261	jfld = jpfld - COUNT( (/.NOT.lhftau/) )
	262	!
	263	! !- allocate the bulk structure
	264	ALLOCATE( sf(jfld), STAT=ierror )
	265	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
	266	DO ifpr= 1, jfld
	267	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
	268	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
[9168]	269	IF( slf_i(ifpr)%nfreqh > 0. .AND. MOD( 3600. * slf_i(ifpr)%nfreqh , REAL(nn_fsbc) * rdt) /= 0. ) &
	270	& CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', &
	271	& ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
[8530]	272
[7163]	273	END DO
	274	! !- fill the bulk structure with namelist informations
	275	CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
	276	!
[7431]	277	IF ( ln_wave ) THEN
[11217]	278	!Activated wave module but neither drag nor stokes drift activated
[9033]	279	IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN
[10425]	280	CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
[11217]	281	!drag coefficient read from wave model definable only with mfs bulk formulae and core
	282	ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR ) THEN
	283	CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
[7431]	284	ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN
[11217]	285	CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
[7431]	286	ENDIF
	287	ELSE
[11217]	288	IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) &
	289	& CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', &
	290	& 'with drag coefficient (ln_cdgw =T) ' , &
	291	& 'or Stokes Drift (ln_sdw=T) ' , &
	292	& 'or ocean stress modification due to waves (ln_tauwoc=T) ', &
	293	& 'or Stokes-Coriolis term (ln_stcori=T)' )
	294	ENDIF
[7431]	295	!
[11217]	296	!
[7163]	297	IF(lwp) THEN !** Control print
	298	!
[11217]	299	WRITE(numout,) ! namelist
[7163]	300	WRITE(numout,*) ' Namelist namsbc_blk (other than data information):'
	301	WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR
	302	WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0
[11284]	303	WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6
[11615]	304	WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF
[7163]	305	WRITE(numout,*) ' add High freq.contribution to the stress module ln_taudif = ', ln_taudif
	306	WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt
	307	WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
	308	WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
	309	WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
	310	WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac
	311	WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
[9019]	312	WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12
	313	WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15
[7163]	314	!
	315	WRITE(numout,*)
	316	SELECT CASE( nblk ) !* Print the choice of bulk algorithm
[9190]	317	CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)'
	318	CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)'
[11284]	319	CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)'
[11615]	320	CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)'
[7163]	321	END SELECT
	322	!
[11111]	323	WRITE(numout,*)
[11615]	324	WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs !LB
	325	WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl !LB
[11284]	326	!
	327	!LB:
[11182]	328	WRITE(numout,*)
[11284]	329	SELECT CASE( nhumi ) !* Print the choice of air humidity
	330	CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]'
	331	CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]'
	332	CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]'
	333	END SELECT
	334	!LB.
[11111]	335	!
[7163]	336	ENDIF
	337	!
	338	END SUBROUTINE sbc_blk_init
	339
	340
[6723]	341	SUBROUTINE sbc_blk( kt )
	342	!!---------------------------------------------------------------------
	343	!! * ROUTINE sbc_blk *
	344	!!
	345	!! ** Purpose : provide at each time step the surface ocean fluxes
[9019]	346	!! (momentum, heat, freshwater and runoff)
[6723]	347	!!
	348	!! ** Method : (1) READ each fluxes in NetCDF files:
	349	!! the 10m wind velocity (i-component) (m/s) at T-point
	350	!! the 10m wind velocity (j-component) (m/s) at T-point
	351	!! the 10m or 2m specific humidity ( % )
	352	!! the solar heat (W/m2)
	353	!! the Long wave (W/m2)
	354	!! the 10m or 2m air temperature (Kelvin)
	355	!! the total precipitation (rain+snow) (Kg/m2/s)
	356	!! the snow (solid prcipitation) (kg/m2/s)
	357	!! the tau diff associated to HF tau (N/m2) at T-point (ln_taudif=T)
	358	!! (2) CALL blk_oce
	359	!!
	360	!! C A U T I O N : never mask the surface stress fields
	361	!! the stress is assumed to be in the (i,j) mesh referential
	362	!!
	363	!! ** Action : defined at each time-step at the air-sea interface
	364	!! - utau, vtau i- and j-component of the wind stress
	365	!! - taum wind stress module at T-point
	366	!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
	367	!! - qns, qsr non-solar and solar heat fluxes
	368	!! - emp upward mass flux (evapo. - precip.)
	369	!! - sfx salt flux due to freezing/melting (non-zero only if ice is present)
	370	!!
	371	!! ** References : Large & Yeager, 2004 / Large & Yeager, 2008
	372	!! Brodeau et al. Ocean Modelling 2010
	373	!!----------------------------------------------------------------------
	374	INTEGER, INTENT(in) :: kt ! ocean time step
	375	!!---------------------------------------------------------------------
	376	!
	377	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
[7163]	378	!
[11266]	379	IF( kt == nit000 ) tsk(:,:) = sst_m(:,:)*tmask(:,:,1) ! no previous estimate of skin temperature => using bulk SST
	380	!
	381	IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) ! compute the surface ocean fluxes using bulk formulea
[6723]	382
	383	#if defined key_cice
	384	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
[7753]	385	qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1)
[11217]	386	IF( ln_dm2dc ) THEN
	387	qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
	388	ELSE
	389	qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
	390	ENDIF
[7753]	391	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)
[11182]	392	!LB:
[11284]	393	SELECT CASE( nhumi )
	394	CASE( np_humi_sph )
	395	qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)
	396	CASE( np_humi_dpt )
	397	qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
	398	CASE( np_humi_rlh )
[11638]	399	qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file
[11284]	400	END SELECT
[11182]	401	!LB.
[7753]	402	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
	403	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
	404	wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
	405	wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
[6723]	406	ENDIF
	407	#endif
	408	!
	409	END SUBROUTINE sbc_blk
	410
	411
	412	SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
	413	!!---------------------------------------------------------------------
	414	!! * ROUTINE blk_oce *
	415	!!
	416	!! ** Purpose : provide the momentum, heat and freshwater fluxes at
	417	!! the ocean surface at each time step
	418	!!
	419	!! ** Method : bulk formulea for the ocean using atmospheric
	420	!! fields read in sbc_read
	421	!!
	422	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
	423	!! - vtau : j-component of the stress at V-point (N/m2)
	424	!! - taum : Wind stress module at T-point (N/m2)
	425	!! - wndm : Wind speed module at T-point (m/s)
	426	!! - qsr : Solar heat flux over the ocean (W/m2)
	427	!! - qns : Non Solar heat flux over the ocean (W/m2)
	428	!! - emp : evaporation minus precipitation (kg/m2/s)
	429	!!
	430	!! ** Nota : sf has to be a dummy argument for AGRIF on NEC
	431	!!---------------------------------------------------------------------
	432	INTEGER , INTENT(in ) :: kt ! time step index
	433	TYPE(fld), INTENT(inout), DIMENSION(:) :: sf ! input data
	434	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pst ! surface temperature [Celcius]
	435	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s]
	436	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s]
	437	!
	438	INTEGER :: ji, jj ! dummy loop indices
	439	REAL(wp) :: zztmp ! local variable
[9019]	440	REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
[11182]	441	REAL(wp), DIMENSION(jpi,jpj) :: zsq ! specific humidity at pst [kg/kg]
[9019]	442	REAL(wp), DIMENSION(jpi,jpj) :: zqlw, zqsb ! long wave and sensible heat fluxes
	443	REAL(wp), DIMENSION(jpi,jpj) :: zqla, zevap ! latent heat fluxes and evaporation
	444	REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin
	445	REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
	446	REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K]
[11182]	447	REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB
[6723]	448	!!---------------------------------------------------------------------
	449	!
[7753]	450	! local scalars ( place there for vector optimisation purposes)
	451	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
[6723]	452
	453	! ----------------------------------------------------------------------------- !
	454	! 0 Wind components and module at T-point relative to the moving ocean !
	455	! ----------------------------------------------------------------------------- !
	456
[7753]	457	! ... components ( U10m - U_oce ) at T-point (unmasked)
[11217]	458	!!gm move zwnd_i (_j) set to zero inside the key_cyclone ???
[7753]	459	zwnd_i(:,:) = 0._wp
	460	zwnd_j(:,:) = 0._wp
[6723]	461	#if defined key_cyclone
	462	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
	463	DO jj = 2, jpjm1
	464	DO ji = fs_2, fs_jpim1 ! vect. opt.
	465	sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
	466	sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
	467	END DO
	468	END DO
	469	#endif
	470	DO jj = 2, jpjm1
	471	DO ji = fs_2, fs_jpim1 ! vect. opt.
	472	zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
	473	zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
	474	END DO
	475	END DO
[10425]	476	CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
[6723]	477	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
[7753]	478	wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
	479	& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
[6723]	480
	481	! ----------------------------------------------------------------------------- !
[11111]	482	! I Solar FLUX !
[6723]	483	! ----------------------------------------------------------------------------- !
	484
	485	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
	486	zztmp = 1. - albo
[11217]	487	IF( ln_dm2dc ) THEN
	488	qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
	489	ELSE
	490	qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)
[6723]	491	ENDIF
	492
	493
	494	! ----------------------------------------------------------------------------- !
	495	! II Turbulent FLUXES !
	496	! ----------------------------------------------------------------------------- !
	497
	498	! ... specific humidity at SST and IST tmask(
[11111]	499	zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
[11182]	500
[11615]	501	IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
[11217]	502	zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action
	503	zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action
	504	END IF
	505
[11182]	506	!LB:
[11284]	507	! zqair = specific humidity of air at "rn_zqt" m above thes sea:
	508	SELECT CASE( nhumi )
	509	CASE( np_humi_sph )
	510	zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity!
	511	CASE( np_humi_dpt )
[11637]	512	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => computing q_air out of d_air and slp !'
[11182]	513	zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
[11284]	514	CASE( np_humi_rlh )
[11637]	515	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => computing q_air out of RH, t_air and slp !'
[11638]	516	zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file
[11284]	517	END SELECT
[11182]	518	!LB.
	519
[6723]	520	!! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
	521	!! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
	522	!! (since reanalysis products provide T at z, not theta !)
[11182]	523	ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt
[6723]	524
[11284]	525
[11615]	526	IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
[11284]	527
[11217]	528	SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point
[11623]	529
[11217]	530	CASE( np_COARE_3p0 )
[11637]	531	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1)
[11617]	532	CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,& ! COARE v3.0
[11615]	533	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
[11266]	534	& Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
[11617]	535	& isecday_utc=nsec_day, plong=glamt(:,:) )
[11623]	536
[11284]	537	CASE( np_COARE_3p6 )
[11637]	538	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => calling "turb_coare3p6" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1)
[11615]	539	CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,& ! COARE v3.6
	540	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
[11284]	541	& Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
[11615]	542	& isecday_utc=nsec_day, plong=glamt(:,:) )
[11284]	543
[11217]	544	CASE( np_ECMWF )
[11637]	545	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1)
[11615]	546	CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & ! ECMWF
	547	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
	548	& Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
[11284]	549
[11217]	550	CASE DEFAULT
[11615]	551	CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' )
[11217]	552	END SELECT
	553
	554	!! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and zsq:
	555	WHERE ( fr_i < 0.001_wp )
	556	! zst and zsq have been updated by cool-skin/warm-layer scheme and we keep it!!!
[11215]	557	zst(:,:) = zst(:,:)*tmask(:,:,1)
	558	zsq(:,:) = zsq(:,:)*tmask(:,:,1)
[11217]	559	ELSEWHERE
	560	! we forget about the update...
	561	zst(:,:) = zqsb(:,:) !LB: using what we backed up before skin-algo
	562	zsq(:,:) = zqla(:,:) !LB: " " "
	563	END WHERE
	564
[11266]	565	!LB: Update of tsk, the "official" array for skin temperature
[11217]	566	tsk(:,:) = zst(:,:)
	567
[11615]	568	ELSE !IF ( ln_skin_cs .OR. ln_skin_wl )
[11217]	569
[11284]	570
[11217]	571	SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point
	572	!
[11615]	573	CASE( np_NCAR )
	574	CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! NCAR-COREv2
	575	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[11217]	576
	577	CASE( np_COARE_3p0 )
[11637]	578	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!'
[11617]	579	CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,& ! COARE v3.0
[11215]	580	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[11111]	581
[11615]	582	CASE( np_COARE_3p6 )
[11637]	583	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!'
[11615]	584	CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,& ! COARE v3.6
	585	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[11217]	586
	587	CASE( np_ECMWF )
[11637]	588	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
[11615]	589	CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & ! ECMWF
[11215]	590	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[11615]	591
[11217]	592	CASE DEFAULT
	593	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
	594	END SELECT
[11111]	595
[11615]	596	END IF ! IF ( ln_skin_cs .OR. ln_skin_wl )
[11111]	597
	598
[11217]	599
	600
[6727]	601	! ! Compute true air density :
[11182]	602	IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove rhoagravrn_zu from SLP...)
	603	rhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) )
[6727]	604	ELSE ! At zt:
[11182]	605	rhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), zqair(:,:), sf(jp_slp)%fnow(:,:,1) )
[6723]	606	END IF
	607
[11217]	608	!! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef.
	609	!! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef.
[7355]	610
[6727]	611	DO jj = 1, jpj ! tau module, i and j component
[6723]	612	DO ji = 1, jpi
[11182]	613	zztmp = rhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed
[6723]	614	taum (ji,jj) = zztmp * wndm (ji,jj)
	615	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
	616	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
	617	END DO
	618	END DO
	619
[6727]	620	! ! add the HF tau contribution to the wind stress module
[7753]	621	IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
[6727]	622
[6723]	623	CALL iom_put( "taum_oce", taum ) ! output wind stress module
	624
	625	! ... utau, vtau at U- and V_points, resp.
	626	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
	627	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
	628	DO jj = 1, jpjm1
	629	DO ji = 1, fs_jpim1
	630	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) &
	631	& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
	632	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) &
	633	& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
	634	END DO
	635	END DO
[10425]	636	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
[6723]	637
	638	! Turbulent fluxes over ocean
	639	! -----------------------------
	640
	641	! zqla used as temporary array, for rho*U (common term of bulk formulae):
[11182]	642	zqla(:,:) = rhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
[6723]	643
	644	IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
	645	!! q_air and t_air are given at 10m (wind reference height)
[11182]	646	zevap(:,:) = rn_efacMAX( 0._wp, zqla(:,:)Ce_atm(:,:)*(zsq(:,:) - zqair(:,:)) ) ! Evaporation, using bulk wind speed
	647	zqsb (:,:) = cp_air(zqair(:,:))zqla(:,:)Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed
[6723]	648	ELSE
	649	!! q_air and t_air are not given at 10m (wind reference height)
	650	! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
[9019]	651	zevap(:,:) = rn_efacMAX( 0._wp, zqla(:,:)Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
[11182]	652	zqsb (:,:) = cp_air(zqair(:,:))zqla(:,:)Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed
[6723]	653	ENDIF
	654
[6727]	655	zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux
[6723]	656
	657
[11111]	658	! ----------------------------------------------------------------------------- !
	659	! III Net longwave radiative FLUX !
	660	! ----------------------------------------------------------------------------- !
	661
[11615]	662	!! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
[11291]	663	zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefanzst(:,:)zst(:,:)zst(:,:)zst(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux
[11111]	664
	665
[6723]	666	IF(ln_ctl) THEN
[9019]	667	CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' )
	668	CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' )
[6723]	669	CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' )
	670	CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' )
	671	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, &
	672	& tab2d_2=vtau , clinfo2= ' vtau : ', mask2=vmask )
	673	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce: wndm : ')
	674	CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce: zst : ')
	675	ENDIF
	676
	677	! ----------------------------------------------------------------------------- !
[11111]	678	! IV Total FLUXES !
[6723]	679	! ----------------------------------------------------------------------------- !
	680	!
[7753]	681	emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.)
	682	& - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1)
	683	!
	684	qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar
[9935]	685	& - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus & ! remove latent melting heat for solid precip
[11111]	686	& - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST !LB??? pst is Celsius !?
[7753]	687	& + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair
	688	& * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp &
	689	& + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
[9935]	690	& * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
[9727]	691	qns(:,:) = qns(:,:) * tmask(:,:,1)
[7753]	692	!
[9570]	693	#if defined key_si3
[9656]	694	qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) ! non solar without emp (only needed by SI3)
[7753]	695	qsr_oce(:,:) = qsr(:,:)
[6723]	696	#endif
	697	!
[11111]	698	!!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN
[11217]	699	CALL iom_put( "rho_air" , rhoa ) ! output air density (kg/m^3) !#LB
	700	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
	701	CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean
	702	CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean
	703	CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean
	704	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
	705	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
	706	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
	707	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
	708	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
	709	CALL iom_put( 'snowpre', sprecip ) ! Snow
	710	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
[11615]	711	IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
[11217]	712	CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius
	713	CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference...
	714	END IF
[11111]	715	!!#LB. ENDIF
[6723]	716	!
	717	IF(ln_ctl) THEN
	718	CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=zqlw , clinfo2=' zqlw : ')
	719	CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ')
	720	CALL prt_ctl(tab2d_1=pst , clinfo1=' blk_oce: pst : ', tab2d_2=emp , clinfo2=' emp : ')
	721	CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, &
	722	& tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask )
	723	ENDIF
	724	!
	725	END SUBROUTINE blk_oce
	726
	727
[9019]	728
[9570]	729	#if defined key_si3
[9019]	730	!!----------------------------------------------------------------------
[9570]	731	!! 'key_si3' SI3 sea-ice model
[9019]	732	!!----------------------------------------------------------------------
	733	!! blk_ice_tau : provide the air-ice stress
	734	!! blk_ice_flx : provide the heat and mass fluxes at air-ice interface
[10534]	735	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
[9019]	736	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
[11217]	737	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
[9019]	738	!!----------------------------------------------------------------------
	739
[6723]	740	SUBROUTINE blk_ice_tau
	741	!!---------------------------------------------------------------------
	742	!! * ROUTINE blk_ice_tau *
	743	!!
	744	!! ** Purpose : provide the surface boundary condition over sea-ice
	745	!!
	746	!! ** Method : compute momentum using bulk formulation
	747	!! formulea, ice variables and read atmospheric fields.
	748	!! NB: ice drag coefficient is assumed to be a constant
	749	!!---------------------------------------------------------------------
	750	INTEGER :: ji, jj ! dummy loop indices
[9019]	751	REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point
	752	REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point
[6723]	753	!!---------------------------------------------------------------------
	754	!
[9019]	755	! set transfer coefficients to default sea-ice values
[11111]	756	Cd_atm(:,:) = rCd_ice
	757	Ch_atm(:,:) = rCd_ice
	758	Ce_atm(:,:) = rCd_ice
[6723]	759
[9019]	760	wndm_ice(:,:) = 0._wp !!gm brutal....
[7355]	761
[9019]	762	! ------------------------------------------------------------ !
	763	! Wind module relative to the moving ice ( U10m - U_ice ) !
	764	! ------------------------------------------------------------ !
[9767]	765	! C-grid ice dynamics : U & V-points (same as ocean)
	766	DO jj = 2, jpjm1
	767	DO ji = fs_2, fs_jpim1 ! vect. opt.
	768	zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )
	769	zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
	770	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
[9019]	771	END DO
[9767]	772	END DO
[10425]	773	CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. )
[9767]	774	!
[9019]	775	! Make ice-atm. drag dependent on ice concentration
	776	IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations
	777	CALL Cdn10_Lupkes2012( Cd_atm )
	778	Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical
	779	ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations
[11217]	780	CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
[7355]	781	ENDIF
	782
[11217]	783	!! CALL iom_put( "rCd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef.
	784	!! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef.
[9019]	785
[6723]	786	! local scalars ( place there for vector optimisation purposes)
	787
[7753]	788	!!gm brutal....
	789	utau_ice (:,:) = 0._wp
	790	vtau_ice (:,:) = 0._wp
	791	!!gm end
[6723]	792
[9019]	793	! ------------------------------------------------------------ !
	794	! Wind stress relative to the moving ice ( U10m - U_ice ) !
	795	! ------------------------------------------------------------ !
[9767]	796	! C-grid ice dynamics : U & V-points (same as ocean)
	797	DO jj = 2, jpjm1
	798	DO ji = fs_2, fs_jpim1 ! vect. opt.
[11182]	799	utau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) &
[9767]	800	& * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
[11182]	801	vtau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) &
[9767]	802	& * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
[6723]	803	END DO
[9767]	804	END DO
[10425]	805	CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
[9019]	806	!
[9767]	807	!
[6723]	808	IF(ln_ctl) THEN
	809	CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ')
	810	CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
	811	ENDIF
[9019]	812	!
[6723]	813	END SUBROUTINE blk_ice_tau
	814
	815
[9019]	816	SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
[6723]	817	!!---------------------------------------------------------------------
	818	!! * ROUTINE blk_ice_flx *
	819	!!
[9019]	820	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
[6723]	821	!!
	822	!! ** Method : compute heat and freshwater exchanged
	823	!! between atmosphere and sea-ice using bulk formulation
	824	!! formulea, ice variables and read atmmospheric fields.
	825	!!
	826	!! caution : the net upward water flux has with mm/day unit
	827	!!---------------------------------------------------------------------
[6727]	828	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature
[9019]	829	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
	830	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
[6727]	831	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
[6723]	832	!!
[6727]	833	INTEGER :: ji, jj, jl ! dummy loop indices
[9454]	834	REAL(wp) :: zst3 ! local variable
[6727]	835	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
[11111]	836	REAL(wp) :: zztmp, z1_rLsub ! - -
[9019]	837	REAL(wp) :: zfr1, zfr2 ! local variables
[9454]	838	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature
[9019]	839	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
	840	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
	841	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
	842	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
[9656]	843	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
[11182]	844	REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB
[6723]	845	!!---------------------------------------------------------------------
	846	!
[11111]	847	zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars
	848	zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
[6723]	849	!
[11182]	850
	851	!LB:
[11284]	852	SELECT CASE( nhumi )
	853	CASE( np_humi_sph )
	854	zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity!
	855	CASE( np_humi_dpt )
[11637]	856	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => ICE !!! computing q_air out of d_air and slp !'
[11182]	857	zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
[11284]	858	CASE( np_humi_rlh )
[11637]	859	IF (lwp) WRITE(numout,) ' ** LB(sbcblk.F90) => ICE !!! computing q_air out of RH, t_air and slp !'
[11638]	860	zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file
[11284]	861	END SELECT
[11182]	862	!LB.
	863
[6723]	864	zztmp = 1. / ( 1. - albo )
[11217]	865	WHERE( ptsu(:,:,:) /= 0._wp )
	866	z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
	867	ELSEWHERE
	868	z1_st(:,:,:) = 0._wp
[9454]	869	END WHERE
[7753]	870	! ! ========================== !
	871	DO jl = 1, jpl ! Loop over ice categories !
	872	! ! ========================== !
[6723]	873	DO jj = 1 , jpj
	874	DO ji = 1, jpi
	875	! ----------------------------!
	876	! I Radiative FLUXES !
	877	! ----------------------------!
[9454]	878	zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
[6723]	879	! Short Wave (sw)
	880	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
	881	! Long Wave (lw)
[11111]	882	z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
[6723]	883	! lw sensitivity
	884	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
	885
	886	! ----------------------------!
	887	! II Turbulent FLUXES !
	888	! ----------------------------!
	889
[9019]	890	! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
[6723]	891	! Sensible Heat
[11182]	892	z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
[6723]	893	! Latent Heat
[11182]	894	qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ch_atm(ji,jj) * wndm_ice(ji,jj) * &
	895	& ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / rhoa(ji,jj) - zqair(ji,jj) ) )
[6723]	896	! Latent heat sensitivity for ice (Dqla/Dt)
	897	IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
[9454]	898	dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) * &
	899	& z1_st(ji,jj,jl)z1_st(ji,jj,jl) EXP(-5897.8 * z1_st(ji,jj,jl))
[6723]	900	ELSE
	901	dqla_ice(ji,jj,jl) = 0._wp
	902	ENDIF
	903
	904	! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
[11182]	905	z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj)
[6723]	906
	907	! ----------------------------!
	908	! III Total FLUXES !
	909	! ----------------------------!
	910	! Downward Non Solar flux
	911	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
	912	! Total non solar heat flux sensitivity for ice
	913	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
	914	END DO
	915	!
	916	END DO
	917	!
	918	END DO
	919	!
[9727]	920	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]
	921	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]
[9019]	922	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
	923	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
[6723]	924
	925	! --- evaporation --- !
[9935]	926	z1_rLsub = 1._wp / rLsub
	927	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation
	928	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT
	929	zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean
[6723]	930
[7753]	931	! --- evaporation minus precipitation --- !
	932	zsnw(:,:) = 0._wp
[9019]	933	CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
	934	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
[7753]	935	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
	936	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
[6723]	937
[7753]	938	! --- heat flux associated with emp --- !
[9019]	939	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
[7753]	940	& + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair
	941	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
[9935]	942	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[7753]	943	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
[9935]	944	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	945
[7753]	946	! --- total solar and non solar fluxes --- !
[9019]	947	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
	948	& + qemp_ice(:,:) + qemp_oce(:,:)
	949	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
[6723]	950
[7753]	951	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
[9935]	952	qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]	953
[7504]	954	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
	955	DO jl = 1, jpl
[9935]	956	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) rcpi * tmask(:,:,1) )
[11217]	957	! ! But we do not have Tice => consider it at 0degC => evap=0
[7504]	958	END DO
	959
[9019]	960	! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
	961	zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! transmission when hi>10cm
	962	zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1
[6723]	963	!
[11217]	964	WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
[9910]	965	qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
[9019]	966	ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm
[9910]	967	qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
[9019]	968	ELSEWHERE ! zero when hs>0
[11217]	969	qtr_ice_top(:,:,:) = 0._wp
[9019]	970	END WHERE
[6723]	971	!
	972	IF(ln_ctl) THEN
	973	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
	974	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
	975	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
	976	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
	977	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
	978	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
	979	ENDIF
	980	!
	981	END SUBROUTINE blk_ice_flx
	982
[11217]	983
[10531]	984	SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
[9019]	985	!!---------------------------------------------------------------------
	986	!! * ROUTINE blk_ice_qcn *
[6723]	987	!!
[9019]	988	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
	989	!! to force sea ice / snow thermodynamics
[10534]	990	!! in the case conduction flux is emulated
[11217]	991	!!
[9019]	992	!! ** Method : compute surface energy balance assuming neglecting heat storage
	993	!! following the 0-layer Semtner (1976) approach
[6727]	994	!!
[9019]	995	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
	996	!! - qcn_ice : surface inner conduction flux (W/m2)
	997	!!
	998	!!---------------------------------------------------------------------
[10531]	999	LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option
[9076]	1000	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
	1001	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
	1002	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
	1003	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
[6723]	1004	!
[9019]	1005	INTEGER , PARAMETER :: nit = 10 ! number of iterations
	1006	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
[6723]	1007	!
[9019]	1008	INTEGER :: ji, jj, jl ! dummy loop indices
	1009	INTEGER :: iter ! local integer
	1010	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
	1011	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
	1012	REAL(wp) :: zqc, zqnet !
	1013	REAL(wp) :: zhe, zqa0 !
	1014	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
	1015	!!---------------------------------------------------------------------
[11217]	1016
[9019]	1017	! -------------------------------------!
	1018	! I Enhanced conduction factor !
	1019	! -------------------------------------!
[10531]	1020	! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
[9019]	1021	! Fichefet and Morales Maqueda, JGR 1997
[6723]	1022	!
[9019]	1023	zgfac(:,:,:) = 1._wp
[11217]	1024
[10531]	1025	IF( ld_virtual_itd ) THEN
[9019]	1026	!
[9935]	1027	zfac = 1._wp / ( rn_cnd_s + rcnd_i )
[9019]	1028	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
	1029	zfac3 = 2._wp / zepsilon
[11217]	1030	!
	1031	DO jl = 1, jpl
[9019]	1032	DO jj = 1 , jpj
	1033	DO ji = 1, jpi
[9935]	1034	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness
[9019]	1035	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
	1036	END DO
	1037	END DO
	1038	END DO
[11217]	1039	!
[10531]	1040	ENDIF
[11217]	1041
[9019]	1042	! -------------------------------------------------------------!
	1043	! II Surface temperature and conduction flux !
	1044	! -------------------------------------------------------------!
[6723]	1045	!
[9935]	1046	zfac = rcnd_i * rn_cnd_s
[6723]	1047	!
[9019]	1048	DO jl = 1, jpl
	1049	DO jj = 1 , jpj
	1050	DO ji = 1, jpi
[11217]	1051	!
[9019]	1052	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
[9935]	1053	& ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
[9019]	1054	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
	1055	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
[9910]	1056	zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
[6727]	1057	!
[9019]	1058	DO iter = 1, nit ! --- Iterative loop
	1059	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
	1060	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
	1061	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
	1062	END DO
	1063	!
	1064	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
	1065	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
	1066	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
[9910]	1067	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) ) &
[11217]	1068	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
[6723]	1069
[9938]	1070	! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
[11217]	1071	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]	1072
[9019]	1073	END DO
[6723]	1074	END DO
[9019]	1075	!
[11217]	1076	END DO
	1077	!
[9019]	1078	END SUBROUTINE blk_ice_qcn
[6723]	1079
[11217]	1080
[7355]	1081	SUBROUTINE Cdn10_Lupkes2012( Cd )
	1082	!!----------------------------------------------------------------------
	1083	!! * ROUTINE Cdn10_Lupkes2012 *
	1084	!!
[11217]	1085	!! ** Purpose : Recompute the neutral air-ice drag referenced at 10m
[9019]	1086	!! to make it dependent on edges at leads, melt ponds and flows.
[7355]	1087	!! After some approximations, this can be resumed to a dependency
	1088	!! on ice concentration.
[11217]	1089	!!
[7355]	1090	!! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50)
	1091	!! with the highest level of approximation: level4, eq.(59)
	1092	!! The generic drag over a cell partly covered by ice can be re-written as follows:
	1093	!!
	1094	!! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)*(nu+1/(10beta)) * A**mu
	1095	!!
	1096	!! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59)
	1097	!! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
	1098	!! A is the concentration of ice minus melt ponds (if any)
	1099	!!
	1100	!! This new drag has a parabolic shape (as a function of A) starting at
[11217]	1101	!! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
[7355]	1102	!! and going down to Cdi(say 1.4e-3) for A=1
	1103	!!
[7507]	1104	!! It is theoretically applicable to all ice conditions (not only MIZ)
[7355]	1105	!! => see Lupkes et al (2013)
	1106	!!
	1107	!! ** References : Lupkes et al. JGR 2012 (theory)
	1108	!! Lupkes et al. GRL 2013 (application to GCM)
	1109	!!
	1110	!!----------------------------------------------------------------------
	1111	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd
	1112	REAL(wp), PARAMETER :: zCe = 2.23e-03_wp
	1113	REAL(wp), PARAMETER :: znu = 1._wp
	1114	REAL(wp), PARAMETER :: zmu = 1._wp
	1115	REAL(wp), PARAMETER :: zbeta = 1._wp
	1116	REAL(wp) :: zcoef
	1117	!!----------------------------------------------------------------------
	1118	zcoef = znu + 1._wp / ( 10._wp * zbeta )
	1119
	1120	! generic drag over a cell partly covered by ice
[7507]	1121	!!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag
[11111]	1122	!! & rCd_ice * at_i_b(:,:) + & ! pure ice drag
[7507]	1123	!! & zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**zmu ! change due to sea-ice morphology
[7355]	1124
	1125	! ice-atm drag
[11111]	1126	Cd(:,:) = rCd_ice + & ! pure ice drag
[7507]	1127	& zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology
[11217]	1128
[7355]	1129	END SUBROUTINE Cdn10_Lupkes2012
[9019]	1130
	1131
	1132	SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
	1133	!!----------------------------------------------------------------------
	1134	!! * ROUTINE Cdn10_Lupkes2015 *
	1135	!!
	1136	!! ** pUrpose : Alternative turbulent transfert coefficients formulation
[11217]	1137	!! between sea-ice and atmosphere with distinct momentum
	1138	!! and heat coefficients depending on sea-ice concentration
[9019]	1139	!! and atmospheric stability (no meltponds effect for now).
[11217]	1140	!!
[9019]	1141	!! ** Method : The parameterization is adapted from Lupkes et al. (2015)
	1142	!! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
	1143	!! it considers specific skin and form drags (Andreas et al. 2010)
[11217]	1144	!! to compute neutral transfert coefficients for both heat and
[9019]	1145	!! momemtum fluxes. Atmospheric stability effect on transfert
	1146	!! coefficient is also taken into account following Louis (1979).
	1147	!!
	1148	!! ** References : Lupkes et al. JGR 2015 (theory)
	1149	!! Lupkes et al. ECHAM6 documentation 2015 (implementation)
	1150	!!
	1151	!!----------------------------------------------------------------------
	1152	!
	1153	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd
	1154	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch
[10511]	1155	REAL(wp), DIMENSION(jpi,jpj) :: ztm_su, zst, zqo_sat, zqi_sat
[9019]	1156	!
	1157	! ECHAM6 constants
	1158	REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m]
	1159	REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m]
	1160	REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m]
	1161	REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41
	1162	REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41
	1163	REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13
	1164	REAL(wp), PARAMETER :: zc2 = zc * zc
	1165	REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14
	1166	REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30
	1167	REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51
	1168	REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56
	1169	REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26
	1170	REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26
	1171	REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma
	1172	REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp
	1173	!
	1174	INTEGER :: ji, jj ! dummy loop indices
	1175	REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu
	1176	REAL(wp) :: zrib_o, zrib_i
	1177	REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice
	1178	REAL(wp) :: zChn_skin_ice, zChn_form_ice
	1179	REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw
	1180	REAL(wp) :: zCdn_form_tmp
	1181	!!----------------------------------------------------------------------
	1182
[10511]	1183	! mean temperature
[11217]	1184	WHERE( at_i_b(:,:) > 1.e-20 )
	1185	ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
	1186	ELSEWHERE
	1187	ztm_su(:,:) = rt0
[10511]	1188	ENDWHERE
[11217]	1189
[9019]	1190	! Momentum Neutral Transfert Coefficients (should be a constant)
	1191	zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40
	1192	zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7
	1193	zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details)
	1194	!zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32)
	1195
	1196	! Heat Neutral Transfert Coefficients
	1197	zChn_skin_ice = vkarmn*2 / ( LOG( rn_zu / z0_ice + 1._wp ) LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) ) ! Eq. 50 + Eq. 52 (cf Lupkes email for details)
[11217]	1198
[9019]	1199	! Atmospheric and Surface Variables
[10511]	1200	zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin
[11111]	1201	zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg]
	1202	zqi_sat(:,:) = q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg] !LB: no 0.98 !!(rdct_qsat_salt)
[9019]	1203	!
	1204	DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility
	1205	DO ji = fs_2, fs_jpim1
	1206	! Virtual potential temperature [K]
[10511]	1207	zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean
	1208	zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice
	1209	zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu
[11217]	1210
[9019]	1211	! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
	1212	zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean
	1213	zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice
[11217]	1214
[9019]	1215	! Momentum and Heat Neutral Transfert Coefficients
	1216	zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40
[11217]	1217	zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53
	1218
[9019]	1219	! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
	1220	z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
	1221	z0i = z0_skin_ice ! over ice (cf Lupkes email for details)
	1222	IF( zrib_o <= 0._wp ) THEN
	1223	zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 10
	1224	zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )*r1_3 / ( Chn_oce(ji,jj) MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26
	1225	& )zgamma )z1_gamma
	1226	ELSE
	1227	zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12
	1228	zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28
	1229	ENDIF
[11217]	1230
[9019]	1231	IF( zrib_i <= 0._wp ) THEN
	1232	zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9
	1233	zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25
	1234	ELSE
	1235	zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11
	1236	zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27
	1237	ENDIF
[11217]	1238
[9019]	1239	! Momentum Transfert Coefficients (Eq. 38)
	1240	Cd(ji,jj) = zCdn_skin_ice * zfmi + &
	1241	& zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
[11217]	1242
[9019]	1243	! Heat Transfert Coefficients (Eq. 49)
	1244	Ch(ji,jj) = zChn_skin_ice * zfhi + &
	1245	& zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
	1246	!
	1247	END DO
	1248	END DO
[10425]	1249	CALL lbc_lnk_multi( 'sbcblk', Cd, 'T', 1., Ch, 'T', 1. )
[9019]	1250	!
	1251	END SUBROUTINE Cdn10_Lupkes2015
	1252
[7355]	1253	#endif
	1254
[6723]	1255	!!======================================================================
	1256	END MODULE sbcblk

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

Download in other formats: