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_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90 @ 11586

Last change on this file since 11586 was 11586, checked in by gsamson, 5 years ago

dev_r11265_ABL : see #2131

merge HPC-13_IRRMANN_BDY_optimization branch @ r11535 (last commit) with dev_r11265_ABL branch @ r11414 (except doc directory)
change ORCA2 results due to changes in HPC-13_IRRMANN_BDY_optimization branch

Property svn:keywords set to Id

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

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

Download in other formats: