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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcblk.F90 in NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90 @ 12588

Last change on this file since 12588 was 12588, checked in by gsamson, 4 years ago

revised ABL model version including:

albmod cleaning
new abl mixing length option (nn_amxl = 3)
use rho_air function from aerobulk everywhere
remove mxl_abl (replaced by master (mxlm_abl) and dissipative (mxlm_abl) mixing lengths)
temporary flag "ln_tpot" to disable potential temperature computation in sbcblk

Property svn:keywords set to Id

File size: 79.4 KB

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

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

Download in other formats: